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ABSTRACT 


The  problem  of  ship  steering  in  canals  and  confined  waters  is  analyzed 
with  emphasis  on  stability  and  bifurcation  analysis.  The  classical  maneuver¬ 
ing  equations  of  motion  augmented  with  a  model  for  ship/canal  interaction 
are  used  to  model  the  open  loop  dynamics.  Coupling  of  a  control  law  and  a 
guidance  scheme  with  appropriate  time  lags  is  employed  to  model  the  essential 
dynamics  of  a  helmsman.  The  complete  system  is  analyzed  using  both  linear 
and  nonlinear  techniques  in  order  to  assess  its  stability  under  finite  distur¬ 
bances.  The  results  indicate  that  for  certain  regions  of  parameters,  limit  cycle 
oscillations  may  develop  which  could  compromise  system  stability  and  safety 
of  operations. 
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I.  INTRODUCTION 


A.  PROBLEM  STATEMENT 

The  problem  of  motion  stability  of  ships  and  other  marine  vehicles  has  been 
the  subject  of  extensive  studies  in  the  past  (Comstock,  1977).  Most  of  these 
studies  are  with  regards  the  ship’s  directional  stability  in  open  waters  under 
open  loop  conditions.  It  is  well  known  that  in  restricted  waters  such  as  canals 
or  rivers,  although  it  is  possible  to  have  positional  stability  in  principle,  in 
reality  it  is  never  the  case.  This  is  due  to  the  destabilizing  effects  of  the  bank 
suction  forces  and  moments.  These  act  always  in  a  destabilizing  fashion;  i.e., 
after  a  small  initial  disturbance  they  produce  a  force  or  a  moment  which  tends 
to  increase  the  ship’s  path  deviation  from  nominal  (Comstock,  1977).  Open 
loop  conditions,  important  as  they  are,  rarely  represent  real  life  applications. 
Ships  traveling  in  canals  are  under  closed  loop  control,  typically  provided 
through  the  helmsman  or  an  autopilot.  It  becomes  then  important  to  assess 
the  stability  of  the  system  under  finite  disturbances  and  incorporating  some 
modeling  of  closed  loop  operations. 


B.  OBJECTIVES  AND  OUTLINE 

This  thesis  utilizes  both  linear  and  nonlinear  techniques  in  order  to  analyze 
the  problem  of  closed  loop  dynamic  stability  of  ships  in  canals.  Ship  modeling 
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is  discussed  in  Chapter  II  and  it  follows  standard  ship  maneuvering  equations 
(Clayton  and  Bishop,  1982)  augmented  by  a  model  for  bank  suction  effects 
(Beck,  1976).  A  Mariner  class  ship  is  selected  as  a  model  because  of  the  avail¬ 
ability  of  data  for  its  hydrodynamic  coefficients  (Comstock,  1977),  (Bernitsas 
and  Kekridis,  1984).  A  heading  control  law  based  on  Nomoto’s  first-order 
model,  coupled  with  a  fine  of  sight  guidance  law  is  chosen  in  order  to  model 
the  fundamental  behavior  of  closed  loop  conditions.  Since,  in  practice  control 
actions  are  hardly  ever  instantaneous,  appropriate  time  lags  are  introduced 
in  the  feedback.  These  are  modeled  using  the  techniques  outlined  in  (Venne, 
1992).  Linear  eigenvalue  analysis  is  performed  in  Chapter  III  in  order  to 
reveal  parametric  stability  boundaries.  It  is  established  that  loss  of  stability 
occurs  when  a  pair  of  complex  conjugate  eigenvalues  crosses  the  imaginary  axis 
which  results  in  periodic  solutions  or  limit  cycles  (Guckenheimer  and  Holmes, 
1983).  A  nonlinear  analysis  based  on  Hopf  bifurcation  methods  (Chow  and 
Mallet-Paret,  1977)  and  (Hassard  and  Wan,  1978)  is  performed  in  Chapter 
IV.  Conclusions  derived  from  this  study  and  some  recommendations  for  fur¬ 
ther  extensions  are  discussed  in  Chapter  V.  The  basic  programs  that  were  used 
to  produce  the  results  are  included  in  the  Appendix. 
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II.  PROBLEM  FORMULATION 


A.  SHIP  DYNAMICS 

Restricting  our  attention  to  the  horizontal  plane,  the  mathenaatical  model 
consists  of  the  nonlinear  sway  (translational  motion  parallel  to  the  vehicle’s 
longitudinal  axis)  and  yaw  (rotational  motion  about  the  vertical  axis  )  equa¬ 
tions  of  motion.  If  we  consider  a  moving  Cartesian  coordinate  frame  located 
at  the  geometric  center  of  the  vehicle,  Newton’s  equations  of  motion  are: 

m{v  +  ur  +  XGr)  =  V,  (1) 

Izr  -|-  mxG^v  +  ur)  —  N  ,  (2) 

where  v  and  r  are  the  relative  sway  and  yaw  velocities  of  the  moving  vehicle 
with  respect  to  the  water,  m  is  the  vehicle’s  mass,  Iz  is  its  moment  of  in¬ 
ertia  with  respect  to  the  vertical  axis,  and  u  is  the  constant  forward  speed. 
Since  all  quantities  will  be  considered  as  dimensionless  in  this  study,  in  the 
following  we  consider  u  to  be  equal  to  one.  xg  is  the  longitudinal  position 
of  the  ship’s  center  of  gravity  with  respect  to  its  centroid,  and  Y,N  represent 
the  total  excitation  sway  force  and  yaw  moment  respectively.  Following  stan¬ 
dard  ship  maneuvering  assumptions,  these  forces  can  be  expressed  as  the  sum 
of  quadratic  drag  terms  and  first-order  memoryless  polynomials  in  v  and  r, 
which  represent  added  mass  and  damping  due  to  water.  In  this  way  they  can 
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be  expressed  by: 


Y  —  YijV  +  YrT  +  YyV  +  -Y^yyV^  +  —YyrrVr^  + 

6  2 

-Yryyvv'^  +  y^r  +  F(V^,  y)  +  YsS  ,  (3) 

N  =  Nyi) Ni-r  +  NyV +'^NyyyV^  ^-]-NyrrVr^  + 

’^N^pyyTv  +  +  A^('0,y)  +  N^S  ,  (4) 

where  Ya,  Na  represent  the  partial  derivatives  with  respect  to  the  indicated 
variable  a  and  6  is  the  effective  rudder  angle.  Y^'tp^y)^  represent  the 

inter action/proximity  forces  and  moments  that  arise  due  to  the  presence  of 
the  canal,  and  shallow  water  effects. 

The  resulting  nonlinear  differential  equations  can  be  non-dimensionalized 
with  respect  to  the  constant  forward  speed  of  the  ship,  u  and  its  length,  1. 
The  dimensionless  time  variable  is  then  equal  to  tu/l.  A  standard  inertial  sys¬ 
tem  (x,  y)  is  introduced  here  where  the  a;— axis  points  in  the  assumed  nominal 
straight  line  path  and  the  y-axis  is  the  distance  from  this  nominal  path,  see 
Figure  1.  Wg  assume  that  the  nominal  straight  line  path  is  along  the  center- 
line  of  the  canal.  In  cases  where  the  concept  of  a  geometric  centerline  is  not 
applicable,  we  assume  that  the  nominal  path  is  along  the  zero  bank  suction 
location.  This  is  consistent  with  recommended  navigation  practices  that  are 
currently  in  use. 

The  complete  model  is  presented  by  the  following  equations: 

(m  -  Yi,)v  -  (Yr  -  mxG)r  =  YyV  +  ^YyyyV^  +  ^YyrrVr^  + 

6  2 
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Figure  1:  Vehicle  geometry  and  definitions  of  symbols 


^rvvTV^  +  {Yr  -  m)T  Y {i> ,  y)  +  ¥56  , 

—  (iVi>  —  mxQ^v  +  {Iz  ““  Nr)T  =  NyV  +  -rNyyyV^  +  -;:‘NyrrVr^+  , 

0  2 

^NryyVV^  +  {Nr  -  mxG)r  +  iV(^,  y)  +  Ns6  , 
ijj  =  r  , 

y  =  sintp  +  V  cos V’  , 


(5) 

(6) 

(7) 

(8) 


where  the  expressions  for  the  ship’s  yaw  rate  as  well  as  its  inertial  position  rate 
have  been  added,  and  ip  is  the  local  heading  angle.  The  interaction/proximity 
forces  and  moments  are  expanded  to  include  up  to  third  order  terms, 


1  1 

Y{ip,y)  =  Y^ip  +  Yyy  +  -Y^^^ip^  +  -Yyyyy^  , 

■^(^)  y)  ~  Nlp'4^  “I"  Nyy  +  ’^N'lpipjlflp  +  '^Nyyyy  , 
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(9) 

(10) 


as  explained  in  more  detail  in  the  folowing  section. 


B.  HYDRODYNAMIC  COEFFICIENTS 

We  chose  a  Mariner  class  ship  as  a  representative  model.  Its  hydrodynamic 
coefficients  and  geometric  properties  were  taken  from  (Comstock,  1977)  and 
(Bernitsas  and  Kekridis,  1984).  Results  from  (Beck,  1976)  were  used  in  order 
to  model  the  bank  suction  forces  and  moments.  Typical  results  are  shown  in 
Figure  2.  These  show  force  and  ftioment  coefficients  versus  ship  deviation  (??), 
and  for  canal  width  w  —  OAL.  Force  and  moment  coefficients  are  nondimen- 
sionalized  with  respect  to  the  water  density  and  the  ship’s  speed  and  length, 
as  is  customary  in  ship  maneuvering.  Since  the  suction  force  and  moment 
must  change  their  sign  as  either  y  oi  ip  changes  its  sign,  they  must  be  modeled 
by  odd  power  pol5Tiomials,  as  was  done  in  the  previous  section.  Numerical 
values  for  the  coefficients  Yy,  Yyyy,  Ny,  and  Nyyy  can  be  found  by  curve  fitting 

of  the  results  of  Figure  2.  Using  a  depth  to  draft  ratio  of  1.9  we  were  able  to 
estimate. 


=  0.02, 

y 

^yyy 

=  0.468, 

Ny 

=  -0.0025, 

N 

■^^yyy 

=  0. 

The  value  of  was  estimated  by  “discretizing”  the  ship  in  two  segments, 
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at  the  bow  and  the  stern,  and  calculating  the  suction  forces  on  each  part 
individually.  Their  resultant  moment  produced, 

=  0.01  . 

The  value  of  was  taken  to  be  zero,  due  to  lack  of  reliable  data.  The  value 

of  was  estimated  to  be  equal  to  — y„  =  0.014  which  is  true  for  motions  along 

the  centerline  of  the  canal  (Comstock,  1977).  Again,  due  to  lack  of  reliable 
data  we  took  =  0. 

C.  CONTROL  LAW 

The  linearized  set  of  equations  (5)  through  (7)  can  be  expressed  in  the 
following  form: 

Ip  =  r  , 

V  —  aix'tp  +  ai2V  +  aizr  +  a^y  +  hi6  ,  (11) 

r  =  a2x'p  +  022V  +  a2zr  +  024?/  +  b28  , 

where  the  coefficients  aij,  bi  are  functions  of  the  vehicle  hydrodynamic  coeffi¬ 
cients  and  geometric  properties  and  are  given  below: 

“11  =  -  Nr)Y^  +  (Yr  -  mxa)N^]  , 

Uyr 

“12  =  -  Nr)Y^  +  {Yr  -  mXG)N^]  , 

Uvr 

“13  =  Ni.){Yr-m)  +  {Yr-rnxG){Nr-rnxG)]  , 

“14  =  -  Ni.)Yy  +  (Yr  -  TnxG)Ny]  , 

Uvr 


021 


=  -  mxG)Y^  +  {m-  Yi)N^]  , 

LJvt 

022  =  “  rnxG)Yy  +  (m  -  yv)iV„]  , 

lyyr 

023  =  -  rnxG)iYr  -m)  +  {m  -  Yy){Nr  -  mxG)]  , 

J-Jvr 

a24  =  -^[{Nv  -  mxG)Yy  +  {m-  Yy)Ny]  , 

Uyr 

bi  =  -^[{Iz-N,)Ys  +  {Yr-mxG)Ns], 

JJvr 

62  =  -  'mxG)Y6  +  {m-  yo)A^5]  , 

J-yyr 

Dw  =  (m  -  Yi,){Iz  -  Nf)  -  (Ni,  -  mxG)iYr  -  mxG)  ■ 

A  control  law  that  could  model  a  human  operator  should  not  be  based  on 
feedback  of  the  side  slip  velocity  v.  Instead,  it  is  more  likely  that  human 
operators  will  respond  to  changes  in  the  ship’s  heading  angle  V'  and  rate  of 
change  of  heading,  r.  Therefore,  we  choose  to  base  our  control  law  on  Nomoto’s 
model,  which  follows  by  assuming  negligible  sway  velocity  u, 

r  =  ar  +  -{-hS  ,  (12) 


where  the  coefficients  are 


O  =  023  , 

b  =  b2, 

c  =  021  . 

A  linear  control  law  can  be  introduced  in  the  form, 

60  =  ki{xl}  -  i^c)  +  hx  ,  (13) 
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where  ipc  is  the  commanded  heading  angle  and  the  gains  fci,  k2  are  computed 
such  that  the  closed  loop  system  (7),  (12),  and  (13)  has  the  desired  dynamics. 
The  existence  of  the  difference  (^  -  V’c)  in  the  control  law  (13)  has  the  effect  of 
trying  to  point  the  ship’s  longitudinal  axis  towards  the  desired  heading.  The 
characteristic  equation  of  the  system  is  obtained  from  (7),  (12),  and  (13)  as 

—  (a  +  hk2)s  —  (c  +  bki)  =  0  . 

If  the  desired  characteristic  equation  is 

5^  +  2<^UnS  +0J^  =  0  , 

the  control  gains  can  be  computed  from 

-t.  =  -^-2. 
b  b 
a  +  2(^ujn 

k2  - 

The  natural  frequency  and  damping  ratio  C  are  selected  based  on  general 
properties  of  second  order  systems.  Relatively  high  values  of  Wn  and  low  values 
of  C  will  correspond  to  a  responsive  operator,  while  the  opposite  is  true  for 
combinations  of  low  and  high  C- 


Finally,  in  order  to  take  into  account  the  effect  of  rudder  saturation,  the 
commanded  rudder  angle  is  given  by 


6  =  5sat  tanh 


(14) 


where  is  the  slope  of  6  at  the  origin  given  by  (13),  and  ^at  is  the  saturation 
limit  on  6,  typically  around  0.4  radians.  The  hyperbolic  tangent  function  is 
used  instead  of  a  hard  saturation  function,  because  of  its  differentiability. 
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D.  GUIDANCE  SCHEME 


Since  the  previous  control  law  stabilizes  the  ship  to  any  commanded  head¬ 
ing  angle,  it  must  be  coupled  with  an  appropriate  orientation  guidance  scheme 
to  provide  path  keeping  along  the  desired  track.  The  simplest  guidance  scheme 
which  models  some  fundamental  aspects  of  helmsman  behavior  is  a  pure  pur¬ 
suit  guidance  where  the  commanded  heading  angle  equals  the  line  of  sight 
angle, 

i^c  =  -  tan"^  ,  (15) 

as  shown  in  Figure  1.  The  ship  is  located  at  {x,y)  and  attempts  to  point  its 
longitudinal  axis  towards  a  target  point  which  is  located  ahead  of  the  ship 
on  the  reference  path  at  a  constant  preview  distance  Xd-  Pursuit  guidance  is 
achieved  by  commanding  a  heading  angle  tpc  equal  to  the  line  of  sight  angle 
(15). 

The  positional  error  information  y  in  equation  (15),  is  assumed  to  lag  the 
actual  error  y  by  an  amount  of  T  seconds,  in  other  words, 

y  =  yit-T).  (16) 

In  this  equation,  the  time  lag  T  models  the  necessary  time  that  it  takes  for 
the  helmsman  to  process  his  path  deviation  and  initiate  appropriate  corrective 
actions. 
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III.  LINEAR  ANALYSIS 


In  this  Chapter,  we  present  a  linearized  analysis  of  the  equations  of  motion. 
The  purpose  of  this  analysis  is  to  assess  stability  or  instability  of  the  equations 
in  response  to  small  deviations  from  the  straight  line  reference  track.  No 
attempt  will  be  made  here  to  characterize  the  transient  response  of  the  system. 
This  can  be  accompHshed  through  numerical  simulations. 


A.  COMBINED  SYSTEM 

If  we  incorporate  the  interaction/proximity  forces  (9)  and  (10)  in  the  ship 
dynamics  model,  equations  (5)  through  (8),  we  end  up  with  the  combined 
system. 

Ip  =  r  , 

(m  —  Yi,)v  -  (Yr  -  mxG)r  =  Y„i;  +  -YwvV^  +  -zYvrrVr^  +  TYrvvfV^-^ 

1  ®  ^ 

(Y-  ~  m)r  +  +  Yyy  +  —Y.^.^'tp^  +  —YyyyU^  +  Ys6  , 

®  1  ®  j 

-{Ni,  -  mxG)v  +  (Iz  -  Nr)r  =  N^v  +  +  -NvrrVr^+  (17) 

1 

-Nrvvfv^  +  {Nr  -  mxG)r  +  N.^ip  +  Nyy  +  -N.^ip^+ 

"^Nyyyy  +  N s6  , 
y  =  sin  Ip  +  V  cos  rp  . 

Study  of  the  asymptotic  properties  of  this  system  is  the  subject  of  this  and 
the  following  chapters. 
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B.  LINEARIZATION 


The  linearized  form  of  equations  (17)  is  the  following, 

Ip  , 

(m  -  Yy)v  -  {Yf  -  mxaY  =  YyV  +  {Y^  -  m)r+ 

—  {Ny  —  mxQ)v  +  [Ig  —  Nr)r  —  NyV  +  (Ny  —  mxG)r+  ' 

iV^V'  +  NyU  +  Ns6  , 
y  =  Ip  +  v  . 

The  rudder  angle  5  has  one  of  the  forms  that  are  developed  below.  The  time 
delay  operator  can  be  expressed  in  terms  of  its  Taylor  series  expansion, 

y{t-T)=y-Ty  +  ^T^y  -  +  . . .  .  (ig) 

Practical  computations  can  be  performed  by  truncating  equation  (19)  to  first, 
second,  or  third  order. 


1-  First  order  approximation  in  y 
We  have, 

y{t-T)^y-Ty  , 


or 


In  its  linear  form. 


and,  therefore. 


y{t  -  T)  —  y  -  T('ip  +  v)  . 


6  =  6o. 


(20) 


Furthermore,  in  its  linear  form  equation  (15)  can  be  written  as 


If  we  incorporate  (20)  into  (13),  we  derive  the  linearized  first  order  approxi¬ 
mation  in  6, 

^  fciT 

6  =  fciV’  +  k^r  H - y - (V’  +  v)  ■  (21) 

Xd  Xd 

2.  Second  order  approximation  in  y 

Keeping  the  second  order  terms  in  y{t  —  T)  we  get, 

y{t-T)=y-Ty  +  ^T^y.  (22) 

If  we  incorporate  (22)  into  (13)  along  with  the  linear  equations  6  —  Sq  and 
=  -y{i  -  T)/Xd,  we  get 

6  =  klip  +  k2r  H — -y - —{'>P  +  v)  +  — (r  -f  v)  .  (23) 

3.  Third  order  approximation  in  y 
In  this  case, 

y{t-T)=y-Ty  +  ,  (24) 

and 

ki  kiT.  ,  kiT^ .  k\T^  . .  ...  ,  . 

8  =  kl'ip  +  k^r  -i - y - {ip  -H  v)  -1-  — — {r  -f-  u)  — - — (r  -I-  v)  .  (25) 

COd  ^d  ^^d  ^^d 

4.  First  order  approximation  in  6 

If  a  time  lag,  T,  exists  on  6  instead  of  simply  y,  then 

6  =  8sa.t  tanh 


where. 


=  6o{t  -T)  =  8o-  T8o  . 
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This  equation  models  a  time  lag  associated  with  the  entire  application  of  the 
necessary  corrective  control  action  and  not  just  its  positional  error.  Using  the 
above  equations  along  with  the  linear  equations  S  =  Si  and  tpc  =  we 

get 

5  =  fciV'  +  +  —y  -  kiTr  -  +  v)  -  k2Tr  .  (26) 

Xd 

5.  First  order  approximation  in  both  y  and  S 

In  a  more  general  setting,  we  can  assume  that  a  time  lag,  Ti,  exists  in  the 
control  law  6,  and  a  different  time  lag,  T2,  is  present  in  the  processing  of  the 
positional  error  y.  Assuming  a  first  order  approximation  for  both,  we  have 

Si  =  6o{t  -  Ti)  =  60-  TiSo  , 


and 

y{t  -  T2)  =  2/  -  T2y  . 


Therefore,  in  this  case  using  the  above  equations  along  with  the  linear  equa¬ 
tions  S  =  61  and  =z  —y(t  —  T2)lxd,  ,we  obtain 


S  —  klip  -\-  k2r  H — -y  —  kiTjr - -^{ip  +v)  — 

Xd  Xd 

kiT2,,  .  kiTiT2, 

- {pp  -b  u)  -| - (r  -b  v)  —  k2Tir  . 

Xd  Xd 


Using  one  of  the  above  expressions,  the  linearized  equations  of  motion  can 


be  written  as 

tp  =  r  , 

V  =  aiiip  -b  ai2V  -b  aisr  -b  ai4y  +  biS  , 

r  =  a2i'lp  +  a22V  -b  023^  -b  (3242/  +  ^>26  , 
y  —  tp  +  v  , 

where  all  coefficients  have  been  previously  defined. 


(28) 
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C.  SERIES  EXPANSIONS  OF  TIME  LAG 


1.  First  order  approximation  in  y 
In  this  case  we  have, 


y(t-T)  =  y  -Ty  . 


Substitution  of  (21)  into  (28)  yields  the 


where, 


1 _ 

1 

- 1 

0 

^21,1 

>^31,1 

1 


0 

^22,1 

^32,1 

1 


■^21,1  —  + 
^22,1  =  0,12  — 
^23,1  =  ®13  + 

^24,1  =  Oi4  + 

■^31,1  =  021  + 

-d.32,1  =  022  — 

^33,1  =  023  + 

■^34,1  =  024  + 


2.  Second  order  approximation  in  y 
In  this  case  we  have. 


y[t-T)=y- 


following  matrix  system. 


1  0  ■ 

■  rp  ' 

^23,1  ■^24,1 

V 

-^33,1  -^34,1 

r 

1 

o 

o 

.  y . 

b\ki  - 


hkiT 

Xd 


hikxT 

Xd 


h\k2 


biki 


Xd 


b2ki  — 


b2kiT 

Xd 


b2kiT 

Xd 


i>2^2 


hki 

Xd 


Ty  +  . 


(29) 
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Substitution  of  (23)  into  (28)  yields  the  following  matrix  system, 


where, 


ij;  ■ 

- 1 

0 

0 

1--1 

0 

1 

■  -0  ■ 

V 

^21,2  ^22,2  ^23,2  ^24,2 

V 

f 

^31,2  ^32,2  ^33,2  ^34,2 

r 

.  y  . 

.  1  1  0  0 

.  y . 

■^21,2 
^22,2 
Ms, 2 

^24,2 

^31,2 

^32,2 

^33,2 

^34,2 


O’liXd  +  bikiXd  —  bikjT 
Xd  -  Q.hbikiT^ 
o-i2^d  —  h\kiT 
Xd  -  0.56ifcir2 
O'lz^d  +  hik^Xd  +  0.5bikiT^ 


Xd  -  0.5bikiT^ 
0'14^d  ■(*  ^1^1 


Xd  -  O.bbikiT^ 


^  b2k^T  ,  62A:iT2  ^ 

Cl2l  T“  ^2^1  I - - - ^21  2 


“22 


Xd  2xd 

b^kiT  ,  b2kiT^ 

-^22,2 


+ 


Xd 


2x„ 


_  ,  ,  b2kiT^  b2kiT^  ^ 

“23  +  b2k2  -i - - 1 - ;; - >123,2 


2xd 


2xri 


„  ,  hki  ,  62^1^2 

“24  H - 1 - :: - ^24  2 


Xd 


2Xa 


3.  Third  order  approximation  in  y 
A  third  order  approximation  in  y  is  based  on. 


(30) 


y(t-T)  =  y-Ty  +  l;T%  -  ir=y<”)  . 

2  6 

Similar  algebraic  steps  as  in  the  previous  two  approximations  result  in  the 
following  eigenvalue  problem. 


0 

0 

0 

0 

1 

01000 

V 

0  0  5333  0  S353 

r 

00010 

y 

-  0  0  B53_3  0  -B55  3 

.  ^2  ^ 
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(31) 


o 

o 

1— 1 

o 

o 

■  V'  ■ 

0  0  0  0  1 

Vi 

^31,3  ^32,3  ^33,3  ^34,3  >^35, 3 

r 

110  0  0 

y 

.  ^51,3  ^52,3  ^53,3  ^54,3  ^55,3  . 

.  V2  . 

where,  =  v,  V2  is  the  rate  of  change  of  v,  and  the  entries  of  the  generalized 
eigenvalue  problem  (31)  are  given  bellow.  Higher  order  approximations  in  T 
can  not  produce  any  usable  stability  results,  since  the  B  matrix  in  (31)  becomes 
singular. 


^31,3  = 

^32,3  = 

^33,3  = 

^34,3  = 

^35,3  = 

■^51 ,3  = 

^52,3  = 

^53,3  = 

>^54, 3  = 

-^55,3  — 

Bzz,z  = 

^35,3  = 

^53, 3  = 


111  +  h\ki  — 
6iA:ir 


hikiT 

Xd 


12  — 


Xd 


13  +  ^>1^2  + 


2xd 


14  + 


hh 

Xd 


-  1 


2xd 
!1  d"  ^2^1  — 

b2kiT 

>2 - 

Xd 

!3  +  ^2^2  + 


b2kiT 

Xd 


bzkiT^^ 

2xd 


+ 


bzki 


•  +  i 
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i 


■655,3 


h2kiT^ 

6xd 


4.  First  order  approximation  in  S 

Substitution  of  (26)  into  (28)  yields  the  following  matrix  system, 
■  1  0  0  O' 

0  1  ^23, 4  0 

0  0  633  4  0 

.00  0  1 _ 

where, 


^21,4 

+  biki 

Xd 

^22,4 

bik{r 

012 

Xd 

^23,4 

= 

013  +  ^1^2  “  b\kiT 

^24,4 

= 

b-^ki 

Ol4  H - 

Xd 

^31,4 

021  +  62^1 

Xd 

CO 

= 

b2kiT 

022 

Xd 

^33,4 

= 

O23  +  ^2^2  ~  i>2^lF 

CO 

= 

62^1 

024  + 

Xd 

623,4 

= 

bik2T 

633,4 

= 

1  +  62^27’ 

5.  First  order  approximation  in  both  y  and  6 

Substitution  of  (27)  into  (28)  yields  the  following  matrix  system, 

■  1  0  0  O' 

0  622,5  623,5  0 

0  632,5  633,5  0 

0  0  0  1 


- 1 

0 

0 

1—* 

0 

_ 1 

■  ^  ■ 

i) 

^21,5  .^22,5  >l23,5  ^24,5 

V 

r 

^31,5  ^32,5  ^33,5  ^34,5 

r 

y  . 

.  1  1  0  0  _ 

.  y 

,  (33) 


^ ' 

1 

o 

T— t 

o 

o 

1 _ 

■  -0  ' 

i) 

•^21,4  ^22,4  ^23,4  >^24,4 

V 

r 

^31,4  ^32,4  ^33,4  ^34,4 

r 

y  . 

.  1  1  0  0  _ 

.  y . 

(32) 
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where, 


•^21,5 

^22,5 

^23,5 

^24,5 

^31,5 

^32,5 

^33,5 

^34,5 

■622,5 

623.5 

632.5 

633.5 


dll  +  biki  — 
bikiTi 


bikiTi  bik\T2 


a\2  — 


Xd 


Xd 

bikiT2 

Xd 


Xd 


®13  +  bik2  —  bikiTi  + 

biki 

«14  H - 

Xd 

tt21  +  ^^2^1  ~ 
b2kiT\ 


bik\TiT2 

Xd 


62^1?!  b2k\T2 


0-22 


Xd 

b2kiT2 


Xd 


Xd 


Xd 


0,23  +  &2^2  —  b2k\Ti  + 

b2k\ 


b2k\T\T2 


024  + 


1  - 


Xd 

b\kiTiT2 


Xd 

b\k2T\ 

b2kiTiT2 

Xd 

1  -l-  62^2^1 


D.  EXACT  COMPUTATION  OF  TIME  LAG 

The  previous  analysis  using  Taylor  series  expansions  of  y{t  —  T)  breaks 
down  for  approximations  beyond  third  order,  as  the  corresponding  generalized 
eigenvalue  problem  becomes  singular.  Thus,  in  order  to  obtain  an  exact  com¬ 
putation  of  the  stability  curves  and  check  the  validity  of  the  calculations,  a 
different  technique  will  be  performed  in  this  section.  This  technique  is  based 
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on  frequency  response  methods  and  it  utilizes  Nyquist’s  criterion  for  stability 
[Friedland  (1986)].We  write  the  system  of  equations  (11)  along  with  equations 
y  =  xp  V  and  6  =  kiip  +  k^r  +  ^y{t  —  T),  in  the  Laplace  domain,  where 

y{t-T) — >ye~^\  (34) 

is  the  time  delay  operator.The  characteristic  equation  of  the  system  is, 

+  ais^  +  0i2S^  +  ass  +  a4  +  +  /?is  +  /?o)— =  0  ,  (35) 

Xd 

where 

ai  =  — ai2  —  023  —  ^>2^2  , 

<^2  =  — oi4  +  O12O23  —  a%thik2  —  022013  +  01262^2  —  ^2^1  —  021  , 

0:3  =  ~®24  —  O13O24  —  0245iA:2  +  O23O14  +  01462^2  —  ^1^1022  — 

O11O22  +  ?>2^lOi2  +  O21O12  , 

a4  =  021O14  +  012024  —  022O14  —  oiia24  +  62^1014  —  hikia2A  , 

P2  =  -biki  , 

Pi  =  —01362^1  +  ^1^1023  —  52^:1  , 

Po  =  012^2^1  ~  ^>1^1022  ~  &2^lOii  +  6lA;i02i  , 


The  characteristic  equation  is  written  as. 


1  +  —G{s)  =  0  , 


where 


,-Ts 


G{s)  =  (^2^  +  P\s  +  /?o)e' 

s'*  +  ais3  +  a2s2  +  ass  +  a4 


(36) 
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is  the  open  loop  transfer  function  ,  and  ^  denotes  the  effective  position  gain. 
For  stability,  we  can  utilize  the  Nyquist  criterion  which  states  that  the  critical 
value  of  Xd  can  be  computed  from, 


—  |(G(ia;)|  =  1  for  axgG{ju)  = -it  (37) 

Xd 


where  j  denotes  the  imaginary  unit.  The  argument,  QxgG(juj)^  of(36)  is, 

,3 


arg  G(joj)  =  —u)T  +  tan 


-1 


_  ,...-1  asuj-aiuj- 
o  /3  2  tan  j  2  .  ' 


The  set  of  equations  (37)  result  in  the  critical  visibility  distance. 


Xd  = 


/3^a;2  +  (/?o  - 


(aso;  —  aiuj^Y  +  “  0^20;^  +  a4)^ 


(39) 


The  value  of  w  in(38)  such  that  argG(juj)  —  — tt  is  the  value  of  the  phase 
crossover  frequency.  After  solving  for  the  phase  crossover  frequency,  the  critical 
value  of  is  obtained  by  direct  substitution  of  this  value  of  lj  into  equation 


(39). 


B.  RESULTS  AND  DISCUSSION 

Typical  results  are  presented  in  Figures  3  through  6.  All  results  shown 
are  nondimensional  unless  otherwise  stated.  The  critical  value  of  Xd  versus  Un 
for  zero  time  lag,  and  parametrized  for  different  values  of  the  damping  ratio 
is  shown  in  Figure  3.  Stability  is  ensured  for  values  of  Xd  greater  than  its 
critical  value.  It  can  be  seen  that  lower  values  of  Un  require  higher  values 
of  Xd  for  stability.  This  means  that  a  less  responsive  helmsman  will  need  a 
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Figure  3:  Critical  xd  versus  con  for  T2  =  0  and  various  values  of  C 


Figure  4:  Critical  xj,  versus  Un  for  C  =  0-8  and  various  values  of  T2 
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longer  lookahead  distance  for  a  stable  operation.  Similar  conclusions  hold  for 
variations  in  the  damping  ratio  (.  In  this  case,  higher  values  of  (  correspond 
to  better  helmsman  response  which  requires  less  lookahead  distance. 

The  effect  of  time  lag  is  shown  in  Figure  4.  Time  lags  are  in  seconds.  It 
can  be  seen  that  reasonable  amounts  of  time  lag  do  not  have  a  serious  effect  on 
stability,  at  least  in  a  linear  sense.  Of  course,  an  amount  of  time  lag  may  have 
a  serious  effect  on  the  transient  response  of  the  system  as  well  as  its  ability  to 
reject  an  external  disturbance.  This  can  only  be  established  by  a  systematic 
series  of  numerical  simulations.  Similar  conclusions  hold  for  non-zero  time 
lags  Ti,  as  shown  in  Figure  5. 

Finally,  the  severe  destabihzing  effect  of  the  canal  is  demonstrated  by  the 
results  of  Figure  6  where  a  comparison  between  canal  and  open  water  results  is 
presented.  It  can  be  seen  that  an  order  of  magnitude  increase  in  the  lookahead 
distance  may  be  required  if  the  same  control  parameters  are  to  be  used  in  both 
open  water  and  a  canal. 
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IV.  NONLINEAR  ANALYSIS 

A.  INTRODUCTION 

It  can  be  numerically  confirmed  that  in  all  cases  of  stability  loss  of  the  pre¬ 
vious  chapter,  one  pair  of  complex  conjugate  eigenvalues  of  the  corresponding 
eigenvalue  problem  crosses  transversally  the  imaginary  axis.  A  situation  like 
this  in  which  a  certain  parameter  is  varied  such  that  the  real  part  of  one  pair 
of  complex  conjugate  eigenvalues  of  the  linearized  system  matrix  crosses  zero, 
results  in  the  system  leaving  its  steady  state  in  an  oscillatory  manner.  This 
loss  of  stability  is  called  Hopf  bifurcation  and  generically  occurs  in  one  of  two 
ways,  supercritical  or  subcritical.  In  the  supercritical  case,  stable  limit  cycles 
are  generated  after  the  nominal  straight  line  motion  loses  its  stability.  The 
amplitudes  of  these  limit  cycles  are  continuously  increasing  as  the  parameter 
distance  from  its  critical  value  is  increased.  For  small  values  of  this  criticality 
distance  the  resulting  limit  cycle  is  of  small  amplitude  and  differs  little  from 
the  initial  nominal  state.  In  the  subcritical  case,  however,  stable  limit  cycles 
are  generated  before  the  nominal  state  loses  its  stability.  Therefore,  depend¬ 
ing  on  the  initial  conditions  it  is  possible  to  diverge  away  from  the  nominal 
straight  line  path  and  converge  towards  a  limit  cycle  even  before  the  nominal 
motion  loses  its  stability.  This  means  that  in  the  subcritical  Hopf  bifurcation 
case  the  domain  of  attraction  of  the  nominal  state  is  decreasing  and  in  fact 
it  shrinks  to  zero  as  the  critical  point  is  approached.  Random  external  dis- 
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turbances  of  sufficient  magnitude  can  throw  the  vehicle  oflF  to  an  oscillatory 
steady  state  even  though  the  nominal  state  may  still  remain  stable.  After  the 
nominal  state  becomes  unstable,  a  discontinuous  increase  in  the  magnitude  of 
motions  is  observed  as  there  exist  no  simple  stable  nearby  attractors  for  the 
vehicle  to  converge  to.  Distinction  between  these  two  qualitatively  different 
types  of  bifurcation  is,  therefore,  essential  in  design.  The  computational  pro¬ 
cedure  requires  higher  order  approximations  in  the  equations  of  motion  and  is 
the  subject  of  the  next  section. 


B.  DETAILED  CALCULATIONS 


The  nonlinear  equations  of  motion  are  written  as. 


i)  —  r  , 

V  =  aiixp  -|-  ai2V  -I-  oisr  -I-  ai4y  +  b-yS'  , 
r  =  a2iV'  +  0.22V  -1-  a23r  -f  a24y  +  b26'  , 
y  =  sin  V’  +  V  cos  V'  , 
where  the  control  law  is, 

6  =  ki'tp  -i-  k2T  4-  ki  tan”’-  ^  , 

Xd 


and,  including  saturation. 


(40) 

(41) 

(42) 

(43) 

(44) 


6  ' 
.  ^sat , 


(45) 


=  ^sat  tanh 

We  perform  a  Taylor  series  expansion  of  the  equations,  keeping  the  first 
vanishing  nonlinear  terms.  Due  to  port/starboard  symmetry  in  the  problem 


non- 
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this  means  expansion  to  third  order  terms, 


(46) 

(47) 


sin  'll;  =  tjj  —  ,  cos  V’  =  1  —  , 


—  1  _  i,/,2 


6' ^6- 


1 


S6l 


sat 


6  =  kiil)  +  k2r  +  —y{t  -T)-  ^y^{t  -  T)  , 

Xd 


(48) 


where  for  consistency,  the  term  y(t  —  T)  is  approximated  by  its  first  order 
expansion  in  T, 


y{t-T)  =  y-Ty  =  y-Ttl;  +  -  Tv  +  .  (49) 

Substitution  of  equations  (46)  to  (49)  into  equations  (40)  to  (45),  produces 
the  system, 

X  =  Ax  +  g(x)  ,  (50) 

where  the  state  variables  vector  is. 


X  =  [il;,v,r,yf  , 


A  is  the  linearized  matrix  at  equilibrium,  and  g(x)  contains  all  third  order 
terms. 

If  T  is  the  matrix  of  eigenvectors  of  A  evaluated  at  the  critical  point  Xd, 
the  linear  change  of  coordinates, 

X  =  Tz  ,  z  =  T“^x  ,  (51) 

transforms  system  (50)  into  its  normal  coordinate  form, 

z  =  T-^ATz  +  T-^g(Tz)  .  (52) 
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At  the  Hopf  bifurcation  point,  matrix  T~^AT  takes  the  form, 

0  —u>Q  0  0  ' 

»ji  — 1  A  rjn  ^0  0  0 

0  0  p  0  ’ 

.0  0  0  g  _ 

where  uq  is  the  imaginary  part  of  the  critical  pair  of  eigenvalues,  and  the 
remaining  two  eigenvalues  p  and  g  are  negative  (or  complex  conjugate  with 
negative  real  parts)  .  Therefore  in  the  new  system  of  coordinates  (51)  the 
dynamics  of  (50)  are  governed  by  a  reduced  two-dimensional  system  zi  and 
Z2  since  the  coordinates  Z3  and  24  correspond  to  the  eigenvalues  p  and  g  and 
are  asymptotically  stable.  For  values  of  close  to  the  bifurcation  point  Xd, 
matrix  T“^AT  is, 

— (wq  0  0 

to  s  (X  s  0  0 

0  0  p  -\-  p'e  0  ’ 

0  0  0  g  -f  g'£  _ 

where,  e  denotes  the  criticality  difference, 

e  =  Xd-  ,  (53) 

to  =  derivative  of  the  real  part  of  the  critical  eigenvalues  with  respect  to  e  , 

a  =  derivative  of  the  imaginary  part  of  the  critical  eigenvalues 

with  respect  to  £  , 

p'  =  derivative  of  p  with  respect  to  e  , 

g'  =  derivative  of  g  with  respect  to  e  . 

Due  to  continuity,  the  eigenvalues  p-t- p'e  and  g-\-g'e  remain  negative  for  small 
nonzero  values  of  e.  Therefore,  the  coordinates  23,  24  correspond  to  negative 


eigenvalues  and  are  asymptotically  stable.  Center  manifold  theory  predicts 
that  the  relationship  between  the  critical  coordinates  zi,  Z2  and  the  stable 
coordinates  Z3,  Z4  is  at  least  of  quadratic  order.  In  fact,  due  to  the  symmetry 
in  our  problem  the  relationship  is  cubic, 

zi  =  0(z3,  z|)  ,  Z2  =  0(zl  z|)  . 

It  follows  that  because  of  this  order,  Z3,  Z4  do  not  influence  the  third  order 
Taylor  series  expansions,  and  can  be  dropped  from  the  equations.  Therefore, 
we  can  write  the  reduced  system  that  describes  the  essential  dynamics  of  (52) 
in  the  form, 

ii  =  a'szi  —  (ct^o  "I"  +  T'i(zi,  Z2)  ,  (54) 

Z2  =  (u;o  +  c>^'e)zi  +  a'ez2  +  F2(zi,Z2)  ,  (55) 

where  Fi,  F2  contain  the  third  order  terms, 


-^l(^l,^2) 

=  Tiizl  +  ri2ZiZ2  +  ri3Ziz^  +  ri4zl  , 

(56) 

F2(zi,  Z2) 

=  r2izf  +  r22zlz2  +  r232lZ2  +  r24Z2  ■ 

(57) 

The  coefficients  are  computable  from  the  previous  Taylor  expansions,  and 
they  are  given  below. 

The  nonlinear  expansion  coefficients  that  are  utilized  in  the  definition 
of  the  cubic  stability  coefficient  /C  are  given  by  the  following; 

rii  =  ni2^2i  +  ni3^3i  +  ni4£4i  +  7112^11  +  7713^21  , 

7’12  =  ni2^22  +  ’^13^32  +  ^14^42  +  ^12^12  +  ’^13^22  j 
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’"13  =  ”12^23  +  ”13^33  +  ’’1443  +  ’ll2<^13  +  ’113(^23  j 

ri4  =  ni2^24  +  ”13^34  +  ”14^44  +  ”12^^14  +  ’ll3C?24  , 

’’21  =  ’’22^21  +  +  ”2441  +  ”22C^11  +  ”23C?21  , 

’’22  =  ”2242  +  ”23^32  +  ”24^2  +  ”22<^12  +  ”23^2  , 

’’23  =  ”22^23  +  ”23^33  +  ”24^43  +  ”22<^13  +  ”23^^23  , 

’’24  =  ”22^4  +  ”23^4  +  ”24^4  +  ”22<^14  +  ”23<^24  ; 

where  we  denote, 

T  =  [rriij]  ,  T~^  =  [n^]  ,  i,j  =  1,...,4. 

For  numerical  purposes  the  critical  eigenvector  of  T  must  be  normalized  so 

that  its  first  nonvanishing  coefficient  is  identically  equal  to  1.  The  coefficients 

£ij  are  given  by, 

~  =  S^ym\im2i  +  6^miim2i  +  4-  5^rr’”ii’”3i 

+^i/)^j/’”ll”l41  +  S.^yymiim^i  +  4t;r’”2l’”31  +  4rr’”2l’”3l 
+4ut/’”2l’”41  +  4yy’”2l’”4i  +  ^rry’”3i’”41  +  4yy’”3l’”4i 
+^^ur’”ll’”2l’”31  +  ^i/'vy’”lim2l’”41  +  5V'ry’”ll’”4l’”31 
”b^t/7'y’”2l’”4l’”31  4“ 

Syvv'^21  4rr’”21  “b  ^yyy’”4i  ) 

•^22  2 

—  =  <5^V’i'(’”'1i’”22  4-  2miimi2m2i)  4-  ^•0uu(’”i2’”2i  +  2miimi2m22) 
4-^^^r(’”ii’”32  +  2miimi2m3i)  4-  <5^rr(’”i2’”3i  4-  2miim3im32) 
4-5^^y(mfim42  4-  2miimi2m4i)  4-  <^t/>yy(’”4i’”i2  +  2miim4im42) 
4-5wr(’”2l’”32  +  2m3im2im22)  4-  4rr(’”22’”3i  4-  2m3im32m2l) 
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+6vvy{'fn\imj^2  +  2m4im2im22)  +  6vyy{m‘22'm\i  +  2m4im42m2i) 

+5rry  (”^3i"^42  +  2m4im3im32)  +  ^ryy(?7l32”^4i  +  2m4im42m3i) 

+5^vr[”^ll^2lW^32  +  rnz\{mnm22  +  mi2m2i)] 

+5^uy[mnm2im42  +  m4i(miim22  +  mi2m2i)] 

+^^ry[”^ll”^4l’^32  +  ’7X3i(miim42  +  ^7112^41)] 

+5„ry[«l2l"l4l"^32  +  m3i(m2im42  +  m22mAl)] 

^i;w(3jT12i^^22)  ”t”  ^rrr(3?7l22772.32)  ^yyy(377l^]^ 77142) 

^  =  ^^^u(^12^21  +  2miimi2m22)  +  +  2mi27772l777,22) 

12^31  +  2miimi2m32)  +  S^r{fniiml2  +  2mi2mz\'mz2) 

(77722^^41  2777x1^^12^^42)  ^^yy (7714277711  “1“  2777x27774177742) 

+5wr(w722”^31  +  2777217712277732)  +  <5vrr(»772lJ7732  +  27773x77722^32) 

+^OTy  (7772277741  +  2777217772277742)  +  ^uyy(7772l77742  +  2777417774277722) 

+^rry  (7773277741  +  2777427773177732)  +  ^ryy  (7773177142  +  2777417774277732) 

[777127772277731  +  77732(777ii77722  +  777i27772l)] 

+^^tjy  [777 127772277741  +  77742  (777ii77722  +  77ll27772l)] 

+c5^ry  [777127774277731  +  77732(777ii77742  +  777i27774i)] 

+^t;j.y  [777227774277731  +  77732(7772177742  +  7772277741)] 

■t■^^^1/)(3777ll?77J^2)  "I”  ^vw(37772i77722)  "H  ^rrr (37773x77732)  "H  ^yyy(37774i77742) 

■7 —  —  777 1277722  ”1”  ^^vtjTTl  1277722  “I”  <5^^j'777]^277732  “f”  ^^rrTTl  12777 32 

Ol 

+6^^^777  ^277742  +  6^yy777  i2777|2  +  6vt;r7772277732  +  6t;rr7772277732 
+6i;^y7772277742  +  6^y7772277742  +  67*7^7773277742  +  6f'yyinz2'^42 
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hi 

b2 

h2 

b2 

hs 

b2 

Ha 

62 

hi 

£42 

£43 

£44 


+6^mi2m22m32  +  6^vy'mi2m22m42  +  S'^rymi2m42ms2  +  6vrym22m42mz2 
+  bjjvv'^TT'22  brrrf^Z2  byyyTn^2  > 


4l 
bi  ’ 
£22 
bi  ’ 
£23 
bi  ’ 
£24 
bi  ’ 


--mil  ^21  +  -mil 


(  ,1  1 

—mil  ^^l2W^2l  +  -miim22  +  -mi2mii 

/  1  1 

— mi2  \  m-iim22  +  -mi2m2i  +  -mi2mii 


--mi2  m22  +  -mi2 


) 

) 


The  coefficients  da  are  given  by, 


*^11  ~  p.  [^ii{Iz  —  Nr)  +  C2i{Yj.  —  mxa)]  , 

UyV 

^12  =  j:^[ci2{Iz  -  Nr)  +  022(54  -  mxo)]  , 
Uyr 

di3  =  -x:  [ciz{Iz  —  Nr)  4-  023(54  —  raxg)]  , 

Uyr 

^14  =  -fj  [ci4{Iz  —  Nr)  +  024(54  —  raxo)]  , 

Uyr 

^21  =  -^[cii(i\4  -  mxo)  +  02i(m  -  54)]  , 

Uyr 

d22  =  7^[ci2(iVi  -  rrixa)  +  022(m  -  Ti,)]  , 

Uyr 

^23  =  -xr-[ciz{Ni  -  mxa)  +  023(m  -  T^,)]  , 

Uyr 

^24  =  7^[ci4(iVi,  -  mxG)  +  024 (m  -  54)]  , 

Uyr 
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where 


Cll 


Cl2 


—  ^^wv'^21  r)  ^’'J'^31^21  "i"  r)^OT^21^31 

0  Z  I 

1 

+ 


+ 


Cl3 


Ci4 


C21 


C22 


C23 


+ 

+ 

+ 


C24 


U 

+  0^yy’^4l  ) 

^y;,„„3m2im22  +  ^vrri'mlim'^  +  2mzimz2m2\) 
]^rvv{iniimz2  +  2m3im2im22)  + 

'^Yyyy3Tn^-jTn42  ) 


+ 


+ 


+ 


X 

+  ■;: 


g^2/yy^^4l^42  j 

^Yvvv^ml2m2i  +  ^Yyrr{ml2m2i  +  2m3im32m22) 

0  z 

-i;.w("r22"^3i  +  2m32m2im22)  +  ^>V’V’V'3"^i2”^n 

o _ 2 _ 


1  2 
D 

g^yy^^42^41  5 

^3^2  ^2 

'qYvvv'^22  2^’'’'^32^22  “I"  ^^^^^22^32 

g^^t/;777-i2  +  Q^yyy'^42  » 

^  3  ^  2  ^  2 
'qNvvv'^21  2 ■^^^^^31^21  2 “^^^^^21^31 

^  AT  3  ^  AT  3 

—  ^'tjj'tprljTnil  +  g-ZVyyym^^  , 


-A/'  ^ 

6 


V’^11  g-^yyy^4i  ? 

..j3m2im22  +  l;Nvrr{mlirn22  +  2mzimz2m2i) 

0  z 

-Nrvv(‘mlimz2  +  2m3im2im22)  +  -N ^^^Zm\-^mi2 

1  2 

-Nyyy3m4:^m42  , 

~iV,,^,u3?7l22^21  "I"  'Z^vrr(jf^Z2^^^  2tHz\1TIZ2'^‘22) 

D  Z 

-Nrvv{‘m22‘mzi  +  2mz2m2im22)  +  -N^MZmX^mn 

2  b 

+  ~Nyyy3Tn  J 


+ 


3  ^  2  ^  2 

vvv'^22  77 32 ^22  “1“  77 ‘^^^‘^^22^ 32 

0  z  z 

1  3  1  3 

'^N'>pil>')(>n^l2  “1“  g'^J/!/2/^42  ■ 
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The  coefficients  in  a  third  order  expansion  of  the  control  law  are  defined 


=  —^{k[)%  +  0.5ki—  +  ^ 

^sat 

C  _  ^  7  /  /7  /  \2  I 

O^ipvv  —  ^2  ^1(^2)  ~  7 

^sat  X^ 

S'iP'ipr  ==  j 

^sat 

^  1  /,  /  vo  r^A)i 


^T/jyy  ““ 


^sat  Xd 

4lt  4  ’ 

^sat 

2  /^2^2  7 


X,  “  ’ 

_  J-t'ti  j.  Iti 

s2  ^^2  2  ■*■  _3  > 

®sat  ^d 

-J-k^h. 

^lat  ^^d  ’ 

--i-A-M 

-^‘^k[k2k2  , 

^sat 

-±2kik',^  -  2^  , 

SL  Xd  xl  ’ 

^sat  Xd 
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and 


***  "  t 


Srrr  — 


'^yyy 


^sat  ' 

1  1 


kiT^ 

3x1 


J_iM 

3  xl 


2ki 

3a:^  ’ 


3x? 


k[  = 


k'o  = 


T 

ki  -  ki—  , 
Xd 

i.  ^ 
-ki—  . 

Xd 


C.  NONLINEAR  COEFFICIENT  K 

If  we  introduce  polar  coordinates  in  the  form, 

zi  =  i? cos 6  ,  Z2  =  RsinO  ,  (58) 

we  can  use  (54)  and  (55)  to  produce  an  equation  describing  the  rate  of  change 
of  the  radial  coordinate  R, 

R  —  oc'eR  +  "PiO^R^  ,  (59) 

and  a  similar  equation  in  the  rate  of  change  of  the  angular  coordinate  9, 

9  =  uj(i  +  u'e  +  Q{9)R?'  .  (60) 

The  system  of  equations  (59)  and  (60)  contains  two  variables,  one  of  which, 
R,  is  slowly  varying  in  time,  whereas  the  other  one,  ^  is  a  fast  variable.  Then, 
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equation  (59)  can  be  averaged  over  one  cycle  in  9  to  produce  an  equation  with 
constant  coefficients  and  similar  stability  properties, 

R  =  a'eR  +  fZR^  ,  (61) 

where, 

^  ^  ^(^)  —  s  (^’’11  ^13  +  ^22  +  3r24)  •  (62) 

Nontrivial  equilibrium  solutions  of  (61)  correspond  to  limit  cycles  in  the 
original  system.  The  nontrivial  equilibrium  of  (61),  Rq,  is  given  by, 

_  I  a'e 

° "  \l~lc  • 

In  our  problem,  since  the  trivial  equilibrium  gains  its  stability  for  Xd  >  a^denticau 
the  coefficient  a'  is  always  negative.  Therefore,  it  can  be  seen  from  (63)  that 
a  limit  cycle  will  exist  provided  K.  and  s  have  the  same  sign.  The  stability 
properties  of  this  limit  cycle  can  be  determined  by  linearization  of  (61)  around 
Rq.  The  linearized  system  is, 

R  =  -2a'e(R  -  Rq)  ,  (64) 

and  its  eigenvalue  is, 

P  =  -2a' e  .  (65) 

We  can  see  that  the  Floquet  exponent  (65)  is  negative  if  e  is  negative,  which 
means  that  K.  must  be  negative.  Therefore,  location  and  stability  of  limit 
cycles  depends  entirely  on  the  cubic  coefficient  1C  which  is  computable  from 
(62) .  We  can  summarize  our  findings  as. 
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Figure  7:  Nonlinear  coefficient  K  versus  Un  for  5sat  =  0.4  and  various  values  of  Q 

•  If  /C  <  0  then  limit  cycles  exist  for  e  <  0  [xd  <  ^^rfcriticai) 
stable. 


•  If  /C  >  0  then  limit  cycles  exist  for  s  >  0  [xd  >  ^dcriticai)  ^^^y 
unstable. 


D.  RESULTS  AND  DISCUSSION 

Typical  results  in  terms  of  the  nonlinear  stability  coefficient  1C  are  shown  in 
Figures  7  through  10.  The  general  conclusion  from  this  series  of  graphs  is  that 
the  bifurcations  are  all  strongly  sub  critical.  This  means  that  any  hnearized 
stability  results  should  be  viewed  with  extreme  caution.  Limit  cycles  exist 
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omega  n 


Figure  8:  Nonlinear  coefficient  K.  versus  a;„  for  C  = 


omega  n 


0.8  and  various  values  of  5sat 


Figure  9;  Nonlinear  coefficient  JC  versus  ujn  for  C 
various  values  of  T2  (sec) 


Figure  10:  Nonlinear  coefficient  K  versus  ujn  with  and  without  channel  effects 

before  stability  in  the  linear  sense  is  lost  and  a  self  sustained  oscillation  in  the 
system  may  develop  as  a  result  of  an  external  disturbance  even  if  the  nominal 
equilibrium  state  is  still  stable.  The  bifurcations  become  more  subcritical; 
i.e.,  K,  is  more  positive,  for  smaller  values  of  Cj  as  Figure  7  shows.  Figure  8 
demonstrates  the  effect  that  the  saturation  limit  5sat  has  on  the  value  of  K. 
Higher  values  of  5sat,  although  are  not  related  to  the  critical  value  of  Xd  result 
in  significant  changes  in  the  value  of  JC.  The  general  trend  is  that  higher  5sat 
results  in  less  subcritical  bifurcations  as  evidenced  by  the  overall  decrease  of 
K.  Figure  9  shows  the  effect  of  non-zero  time  lag  on  the  nonlinear  stability 
coefficient.  It  can  be  seen  that,  like  the  linear  stability  results,  the  effect  is 
Tninimal.  Finally,  Figure  10  shows  the  canal  effect  on  K,.  This  figure  was 
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produced  for  zero  time  lag,  C  =  0.8,  and  Ssat  =  0.4.  It  can  be  seen  that  the 
existence  of  a  canal  causes  the  bifurcations  to  be  much  more  sub  critical  than 
the  open  water  case.  This  demonstrates  the  severe  destabilizing  effect  that  the 
canal  introduces  in  both  the  linearized  and  the  nonlinear  analysis. 


42 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  thesis  presented  a  comprehensive  study  of  linear  and  nonlinear  stability 
properties  of  straightline  motion  of  surface  ships  in  confined  waters.  The 
classical  maneuvering  equations  of  motion  incorporating  canal  suction  effects, 
were  coupled  with  appropriate  navigation,  gauidance,  and  control  laws  in  order 
to  mimic  the  helmsman’s  behavior.  The  main  conclusions  from  this  study  can 
be  summarized  as  follows: 

1.  There  exists  a  critical  preview  distance  for  straightline  positional  stabil¬ 
ity.  For  values  less  than  the  critical  distance,  the  system  is  unstable. 

2.  Including  canal  effects,  the  critical  preview  distance  may  be  an  order  of 
magnitude  higher  than  in  open  water.  If  the  same  preview  distance  is  to 
be  used  in  both  cases,  the  corresponding  control  law  for  canal  maneuver¬ 
ing  must  be  considerably  more  responsive  than  in  open  waters. 

3.  The  critical  preview  distance  is  monotonically  decreasing  for  increasing 
control  law  responsiveness.  This  means  that  in  order  to  accomodate 
smaller  values  of  the  preview  distance,  more  responsive  control  laws  are 
required  . 

4.  Physically  realizable  time  lags  do  not  seem  to  have  a  significant  effect  on 
the  values  of  the  preview  distance. 

5.  As  the  preview  distance  becomes  less  than  its  critical  value,  one  pair  of 
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complex  conjugate  eigenvalues  of  the  linearized  system  matrix  crosses 
the  imaginary  axis.  This  corresponds  to  a  bifurcation  to  periodic  so¬ 
lutions  (Hopf  bifurcation)  and  the  system  exhibits  oscillatory  behavior, 
also  known  as  limit  cycles. 

6.  Higher  order  approximations  in  the  equations  of  motion  were  utilized  in 
order  to  assess  the  stability  of  the  resulting  limit  cycles.  It  was  found 
that  in  all  cases,  the  limit  cycles  were  unstable.  This  has  the  following 
implications: 

(a)  It  is  possible  for  the  system  to  lose  its  stability  even  before  the  critical 
preview  distance  is  crossed.  This  means  that  all  linearized  stability 
analysis  results  should  be  viewed  with  extreme  caution. 

(b)  As  the  critical  preview  distance  is  crossed,  it  is  expected  that  the 
system  will  develop  limit  cycles  of  large  amplitude.  This  clearly 
presents  a  dangerous  situation  which  should  be  avoided  in  practice, 
by  appropriate  changes  in  the  design  parameters. 

Recommendations  for  further  research  include  the  following: 

1.  Simulation  studies  in  order  to  verify  the  limit  cycle  behavior. 

2.  Study  of  different  ship  characteristics  and  canal  geometry. 
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APPENDIX 


The  following  is  a  list  and  description  of  the  computer  programs  used  in  this 
thesis.  The  programs  are  written  in  FORTRAN.  Complete  printouts  of  the 
programs  follow.  Standard  eigenvalue/eigenvector  numerical  analysis  subrou¬ 
tines  are  required  for  all  programs. 

•  THESISl.FOR 

Calculation  of  critical  Xd-  First  order  approximation  for  time  lag  T2. 

•  THESIS2.FOR 

Calculation  of  critical  Xd-  Second  order  approximation  for  time  lag  T2. 

•  THESIS3.FOR 

Calculation  of  critical  Xd-  Third  order  approximation  for  time  lag  r2. 

•  THESIS4.FOR 

Calculation  of  critical  Xd-  First  order  approximation  for  time  lag  Ti. 

•  THESIS5.FOR 

Calculation  of  critical  Xd-  First  order  approximation  for  time  lags  Ti  and 

T2. 

•  HOPF.FOR 

Calculation  of  the  nonlinear  cubic  stabihty  coefficient  K,.  Requires  the 

output  of  one  of  the  previous  programs  as  its  input. 
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C  PROGRAM  THESISl.FOR  (Time  Delay-lst  Order  Approx.  T2) 

C  BIFURCATION  ANALYSIS 

C 

PROGRAM  THESIS 1 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2.L,NR,NV.NDELTA,NPSI,NY,IZ,MASS, 
&  NRDOT.NVDOT 

DIMENSION  A(4,4) .FV1(4) ,IV1(4) ,ZZZ(4,4) ,WR(4) .WI(4) 

C 

OPEN  (11,FILE=’BIF1.RES’) 

OPEN  (12,FILE=>BIF2.RES’) 

OPEN  (13,FILE=’BIF3.RES’) 
open  (15,file=’eig.resl’) 

C 


Vehicle  Paraineters 

IZ 

=0.0 

L 

=528 

RHO 

=1.94 

XG 

=0.0 

MASS 

=0.0088 

U 

ft 

o 

YRDOT 

=  0.00000 

YVDOT 

=-0.00912 

YR 

=+0,00456 

YV 

=-0.01434 

YPSI 

=  0.01400 

YY 

=  0.02000 

YDELTA=  0.00278 
NRDOT  =-0.00115 
NVDOT  =  0.00000 
NR  =-0.00296 
NV  =-0.00460 
NPSI  =  0.01000 
NY  =-0.00250 
NDELTA=-0 . 00139 
WRITE  (*.1001) 

READ  (*,*)  WNMIN.WNMAX.IWN 
WRITE  (*,1002) 

READ  (*,*)  XDMIN.XDMAX.IXD 
WRITE  (*.1003) 

READ  (*,*)  ZETA 
WRITE (*,1100) 

READ  (*,*)  TL 
TL=TL*U/L 
C 

DVR  =(IZ-NRDOT)*(MASS-YVDOT)- 
&  (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 
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AA11=((IZ-NRD0T) *YPSI-(MASS*XG-YRDOT)*NPSI)/DVR 
AA12=( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV)/DVR 
AA21=( (NVDOT-MASS*XG) *YPSI+(MASS-YVDOT) *NPSI) /DVR 
AA22=( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13= ( ( IZ-NRDOT) * ( YR-MASS) + 

&  (YRDOT-MASS*XG) ♦ (NR-MASS*XG) ) /DVR 

AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

&  (MASS-YVDOT) *(NR-MASS*XG)) /DVR 

AA14=( (IZ-NRDOT) *YY+(YRDOT-MASS*XG) *NY)/DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 
BBl  = ( (IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) ♦YDELTA) /DVR 

AN0M=AA23 
BN0M=BB2 
CN0M=AA21 
EPS  =l.D-5 
ILMAX=1500 

DO  1  1=1, IWN 

WRITE  (*,2001)  I, IWN 

WN=WNMIN+ ( I-l) * (WNMAX-WNMIN) / (IWN- 1) 

Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- ( ANOM+2 . DO*ZETA*WN) /BNOM 
DO  2  J=1,IXD 

XD=XDMIN+( J-1) * (XDMAX-XDMIN) / (IXD-1) 

CALL  LINEAR(K1 , K2 , XD , AAll , AA12 , AA13 , AA14 , AA21 , AA22 , AA23 , 
&  AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 
write(15,*)DE0S,XD,WN 

IF  (J.GT.l)  GO  TO  10 

DEOSOO=DEOS 

XDOO  =XD 

LL=0 

GO  TO  2 

10  DEOSNN=DEOS 
XDNN  =XD 
PR=DEOSNN*DEOSOO 
IF  (PR.GT.O.DO)  GO  TO  3 
LL=LL+1 

IF  (LL.GT.3)  STOP  1000 
IL=0 

XDO=XDOO 

XDN=XDNN 

DEOSO=DEOSOO 
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DEOSN=DEOSNN 
6  XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 
C 

CALL  LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA2 1 , 
&  AA22 , AA23 , AA24 , BBl , BB2 , A , TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL (DECS , WR , WI , FREQ) 

DEOSM=DEOS 


XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS (XDL-XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 


GO  TO  4 

5  IF  (PRR.GT.O.DO)  STOP  3200 
XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 

4  LLL=10+LL 

WRITE  (LLL,*)  XD,WN 
3  XDOO=XDNN 

DEOSOO=DEOSNN 
2  CONTINUE 
1  CONTINUE 


C 


1001  FORMAT 

1002  FORMAT 

1003  FORMAT 
1100  FORMAT 
2001  FORMAT 


('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN 
(^  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD 
ENTER  DAMPING  RATIO’) 

(’  ENTER  TIME  LAG  TL  (dimensional)’) 
(215) 


(dimensionless) ’) 
(dimensionless) ’) 
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END 


C 

SUBROUTINE  DSTABL (DECS, WR.WI, OMEGA) 

C 

C  EVALUATES  THE  EIGENVALUE  WITH  THE  MAXIMUM  REAL  PART 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  WR(4).WI(4) 

DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I) .LT.DEOS)  GO  TO  1 
DEOS=WR(I) 

IJ=I 

1  CONTINUE 
OMEGA=WI(IJ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 

C 

SUBROUTINE  LINEAR (K1 , K2 , XD , AAll , AA12 , AA13 , AA14 , AA2 1 , 

&  AA22,AA23,AA24,BB1,BB2,A,TL) 

C 

C  FORMS  THE  LINEARIZED  MATRIX  A  (time  delay  1st  order  approximation 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2 
DIMENSION  A (4, 4) 

A(l.l)=O.ODO 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(2,1)=AA11+BB1*K1-BB1*K1*TL/XD 

A(2 , 2) =AA12-BB1*K1*TL/XD 

A(2,3)=AA13+BB1*K2 

A(2,4)=AA14+BB1*K1/XD 

A (3 . 1) =AA21+BB2*K1-BB2*K1*TL/XD 

A (3 , 2) =AA22-BB2*K1*TL/XD 

A(3,3)=AA23+BB2*K2 

A (3 , 4) =AA24+BB2*K1/XD 

A(4,1)=1.0DO 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

RETURN 

END 
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C  PROGRAM  THESIS2.F0R  (Time  Delay- 2iid  Order  Approx  T2) 

C  BIFURCATION  ANALYSIS 
C 

PROGRAM  THESIS2 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2,L,NR.NV,NDELTA,NPSI.NY,IZ.MASS. 
&  NRDOT.NVDOT 

DIMENSION  A(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) ,WR(4) ,WI(4) 

C 

OPEN  (11,FILE=’BIF1.RES’) 

OPEN  (12,FILE='BIF2.RES>) 

OPEN  (13,FILE=>BIF3.RES’) 

open  (15,file=’eig.resl’) 

C 

C  Vehicle  Parameters: 

IZ  =0.0 
L  =528 
RHO  =1.94 
XG  =0.0 
MASS  =0.0088 
U  =1.0 
C 

YRDOT  =  0.00000 
YVDOT  =-0.00912 
YR  =+0.00456 
YV  =-0.01434 

YPSI  =  0.01400 
YY  =  0.02000 
YDELTA=  0.00278 

NRDOT  =-0.00115 
NVDOT  =  0- 00000 
NR  =-0.00296 
NV  =-0.00460 

NPSI  =  0.01000 
NY  =-0 . 00250 

NDELTA=-0. 00139 
WRITE  (*,1001) 

READ  (*,*)  WNMIN,WNMAX,IWN 
WRITE  (*,1002) 

READ  ( * , * )  XDMIN , XDMAX , IXD 
WRITE  (*,1003) 
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READ  (*,*)  ZETA 
WRITE(*,1100) 

READ  (*,*)  TL 
TL=TL*U/L 
C 

DVR  =(IZ-NRDOT)*(MASS-YVDOT)- 
&  (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

AA11= { (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12= ( (IZ-NRDDT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21=( (NVDOT-MASS*XG) *YPSI+(MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13= ( (IZ-NRDOT) * ( YR-MASS) + 

&  (YRDOT-MASS*XG)*(NR-MASS*XG))/DVR 

AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

&  (MASS-YVDOT) *(NR-MASS*XG)) /DVR 

AA14= ( ( IZ-NRDOT) *YY+ (YRDOT-MASS+XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BBl  = ( ( IZ-NRDOT) ♦YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2  =( (MASS-YVDOT) ♦NDELTA-(MASS*XG-NVDOT)*YDELTA) /DVR 
C 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS  =l.D-5 
ILMAX=1500 
C 

DO  1  I=l,IVfN 

WRITE  (*,2001)  I.IWN 
WN=WNMIN+(I-1) * (WNMAX-WNMIN) / (IWN-1) 

Kl=- ( WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- ( ANOM+2 . DO*ZETA*WN) /BNOM 

DO  2  J=1.IXD 

XD=XDMIN+ ( J- 1 ) * (XDMAX-XDMIN) / ( IXD- 1 ) 

C 

CALL  LINEAR(K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA21 . AA22 , AA23 , 
&  AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WB,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 
write ( 15 , *) DEOS , XD , WN 
C 

IF  (J.GT.l)  GO  TO  10 

DEOSOO=DEOS 

XDOO  =XD 
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LL=0 
GO  TO  2 

10  DEOSNN=DEOS 
XDNN  =XD 
PR=DEOSNN*DEOSDO 
IF  (PR.GT.O.DO)  GO  TO  3 
LL=LL+1 

IF  (LL.GT.3)  STOP  1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
6  XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL  LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 ,  AA2 1 , 
k  AA22,AA23,AA24,BB1,BB2,A,TL) 

CALL  RG (4 , 4 , A , WR , WI , 0 , ZZZ . I VI , FVl , lERR) 

CALL  DSTABL (DECS, WR.WI, FREQ) 

DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDL-XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 

GO  TO  4 

5  IF  (PRR.GT.O.DO)  STOP  3200 
XDD=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS (XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 
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4  LLL=10+LL 

WRITE  (LLL,*)  XD,WN 
3  XDOO=XDNN 

DEOSOO=DEOSNN 
2  CONTINUE 
1  CONTINUE 


C 


1001 

1002 

1003 

1100 

2001 


FORMAT  (’  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN 
FORMAT  (’  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD 
FORMAT  ('  ENTER  DAMPING  RATIO’) 

FORMAT  (’  ENTER  TIME  LAG  TL  (dimensional)’) 
FORMAT  (215) 

END 


(dimensionless)  ’) 
(dimensionless)  ’) 


C 

SUBROUTINE  DSTABL(DEOS,WR,WI , OMEGA) 

C 

C  EVALUATES  THE  EIGENVALUE  WITH  THE  MAXIMUM  REAL  PART 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  WR(4),WI(4) 

DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I) .LT.DEOS)  GO  TO  1 
DEOS=WR(I) 

IJ=I 

1  CONTINUE 
OMEGA=WI(IJ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 

C 

SUBROUTINE  LINEAR (K1 , K2 , XD , AAll , AA12 , AA13 , AA14 , AA2 1 , 

&  AA22 , AA23 , AA24 , BBl , BB2 , A , TL) 

C 

C  FORMS  THE  LINEARIZED  MATRIX  A  (time  delay  1st  order  approximation) 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2 
DIMENSION  A (4, 4) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=O.ODO 

A(2 , 1) = (AA11*XD+BB1*K1*XD-BB1*K1*TL) / (XD-0 . 50D0*BB1*K1*TL*TL) 

A (2 , 2) = (AA12*XD-BB1*K1*TL) / (XD-0 . 50DO*BB1*K1*TL*TL) 

A (2 ,3) = (AA13*XD+BBl*K2*XD+0 . 50D0*BB1*K1*TL*TL) / 

&  (XD-0 . 50D0*BB1*K1*TL*TL) 

A (2 , 4) = (AA14*XD+BB1*K1) / (XD-0 . 50D0*BB1*K1*TL*TL) 
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A(3,1)=AA21+BB2*K1-BB2*K1*TL/XD+(BB2*K1*TL*TL/(2.0D0*XD))*A(2,1) 

A (3 , 2) =AA22-BB2*K1*TL/XD+ (BB2*K1*TL*TL/ (2 . ODO*XD) ) *A (2 , 2) 

A (3 , 3) =AA23+BB2*K2+ (BB2*K1*TL*TL/ (2 . ODO*XD) ) + 

&  (BB2*K1*TL*TL/(2.0*XD))*  A (2, 3) 

A (3 , 4) =AA24+BB2*K1/XD+ (BB2*K1 *TL*TL/ (2 . ODO*XD) ) * A (2 , 4) 

A(4,1)=1.0D0 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

RETURN 

END 
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C  PROGRAM  THESIS3.F0R  (Time  Delay-3rd  Order  Approx  T2) 

C  BIFURCATION  ANALYSIS 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2,L,NR,NV,IZ,MASS,NDELTA.NPSI ,NY, 
&  NRDOT.NVDOT 

DIMENSION  A(5,5) ,B(5,5) ,FV1(5) ,IV1(5) ,ZZZ(5,5) ,ALFR(5) , 
&  ALFI(5),BETA(5),WR(5),WI(5) 

C 

OPEN  (11,FILE=’BIF1.RES’) 

OPEN  (12,FILE=’BIF2.RES’) 

OPEN  (13,FILE=’BIF3.RES’) 


Vehicle  Parameters 

IZ 

=0.0 

L 

=528 

RHO 

=1.94 

XG 

=0.0 

MASS 

=0.0088 

U 

=1.0 

YRDOT 

=  0.00000 

YVDOT 

=-0.00912 

YR 

=+0.00456 

YV 

=-0.01434 

YPSI 

=  0.01400 

YY 

=  0.02000 

yDELTA=  0.00278 
NRDOT  =-0.00115 
NVDOT  =  0.00000 
NR  =-0.00296 
NV  =-0.00460 
NPSI  =  0.01000 
NY  =-0.00250 
NDELTA=-0.00139 
WRITE  (♦,1001) 

READ  (♦,♦)  WNMIN,WNMAX,IWN 
WRITE  (♦,1002) 

READ  (♦,*)  XDMIN,XDMAX,IXD 
WRITE  (^,1003) 

READ  (♦,♦)  ZETA 
WRITE(^,1100) 

READ  (♦,♦)  TL 
TL=TL^U/L 
C 

DVR  =(IZ-NRDOT)^(MASS-YVDOT)- 
&  (MASS^XG-YRDOT)^(MASS^XG-NVDOT) 
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AA11= ( (IZ-NRDOT) *YPSI- (MASS+XG-YRDOT) *NPSI) /DVR 
AA12= ( (IZ-NRDOT) *YV- (MASS+XG-YRDOT) *NV) /DVR 
AA21= ( (NVDDT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13= ( ( IZ-NRDOT) * ( YR-MASS) + 

&  (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 

AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

&  (MASS-YVDOT) ♦(NR-MASS*XG)) /DVR 

AA14= ( (IZ-NRDOT) ♦YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BBl  = ( (IZ-NRDOT) ♦YDELTA- (MASS*XG-YRDOT) ♦NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS  =l.D-5 
ILMAX=1500 
C 

DO  1  1=1, IWN 

WRITE  (*,2001)  I, IWN 

WN=WNMIN+ ( I- 1 ) * (WNMAX-WNMIN) / ( IWN- 1 ) 

Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- (ANOM+2 . DO*ZETA*WN) /BNOM 
DO  2  J=1,IXD 

XD=XDMIN+ ( J- 1 ) * (XDMAX-XDMIN) / ( IXD- 1 ) 

C 

CALL  LINEAR (K1 , K2 , XD , AAl 1 , AAl 2 , AA13 , AA14 , 

&  AA21,AA22,AA23,AA24,BB1,BB2,A,B,TL) 

CALL  RGG (5 , 5 , A . B . ALFR , ALFI , BETA, 0 , ZZZ , lERR) 

DO  11  IJE=1,5 

WR(IJE)=ALFR(IJE)/BETA(IJE) 

WI ( IJE) =ALFI ( I JE) /BETA ( I JE) 

11  CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

C 

IF  (J.GT.l)  GO  TO  10 

DEDSOD=DEOS 

XDOO  =XD 

LL=0 

GO  TO  2 

10  DEDSNN=DEOS 

XDNN  =XD 
PR=DEOSNN*DEOSOO 
IF  (PR.GT.O.DO)  GO  TO  3 
LL=LL+1 
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IF  (LL.GT.3)  STOP  1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEDSN=DEOSNN 
6  XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL  LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA21 , AA22 , 
&  AA23 , AA24 , BB 1 , BB2 , A , B , TL) 

CALL  RGG (5 , 5 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , lERR) 

DO  12  IJE=1,5 

WR ( I JE) -ALFR ( I JE) /BETA ( IJE) 

WI ( I JE) =ALFI ( I JE) /BETA ( I JE) 

12  CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

DEOSM-DEOS 

XDM-XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

XDO-XDL 

XDN-XDM 

DEOSO-DEOSL 

DEOSN-DEOSM 

IL-IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF-DABS(XDL-XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD-XDM 

GO  TO  4 

5  IF  (PRR.GT.O.DO)  STOP  3200 

XDO-XDM 
XDN-XDR 
DEOSO-DEOSM 
DEOSN-DEOSR 
IL-IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO  6 
XD-XDM 

4  LLL-IO+LL 

WRITE  (LLL,*)  XD,WN 
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3  XDOO=XDNN 

DEOSOO=DEOSNN 
2  CONTINUE 
1  CONTINUE 
C 

1001  FORMAT  (»  ENTER  MIN.  MAX,  AND  INCREMENTS  OF  WN  (dimensionless)') 

1002  FORMAT  (’  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD  (dimensionless)') 

1003  FORMAT  ('  ENTER  DAMPING  RATIO') 

1100  FORMAT  ('  ENTER  TIME  LAG  TL  (dimensional)') 

2001  FORMAT  (215) 

END 

C 

SUBROUTINE  DSTABL (DEOS , WR , WI , OMEGA) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  WR(5),WI(5) 

DEOS=-1.0D+20 
DO  1  1=1,5 

IF  (WR(I) .LT.DEOS)  GO  TO  1 
DEOS=WR(I) 

IJ=I 

1  CONTINUE 
OMEGA=WI(IJ) 

OMEGA=D ABS ( OMEGA ) 

RETURN 

END 

C 

SUBROUTINE  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23,AA24 
&  ,BB1,BB2,A,B,TL) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2 
DIMENSION  A(5,5),B(5,5) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=O.ODO 

A(1,5)=0.0D0 

A(2,1)=O.ODO 

A(2,2)=O.ODO 

A(2,3)=0.0D0 

A(2,4)=0.0D0 

A(2,5)=1.0D0 

A(3, 1)=AA11- (BB1*K1*TL/XD) +BB1*K1 
A (3 , 2) =AA12-BB1*K1*TL/XD 

A (3 , 3) =AA13+BB1*K2+ (BB1*K1*TL*TL/ (2 . DO*XD) ) 

A(3,4)=AA14+BB1*K1/XD 
A (3 , 5) =-l .D0+ (BB1*K1*TL*TL/ (2 . DO*XD) ) 

A(4,1)=1.0D0 
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A(4,2)=1.0DO 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

A(4,5)=0.0D0 

A(5, 1)=AA21+BB2*K1-BB2*K1*TL/XD 
A (5 . 2) =AA22-BB2*K1*TL/XD 

A (5 , 3) =AA23+BB2*K2+ (BB2*K1 *TL*TL/ (2 . DO*XD) ) 

A (5 , 4) =AA24+BB2*K1/XD 

A(5 , 5)= (BB2*K1*TL*TL/ (2 .DO*XD) ) 

C 

B(1,1)=1.0D0 

B(1,2)=0.0D0 

B(1,3)=O.ODO 

B(1,4)=O.ODO 

B(1,5)=0.0D0 

B(2,1)=O.ODO 

B(2,2)=1.0DO 

B(2,3)=O.ODO 

B(2,4)=O.ODO 

B(2,5)=O.ODO 

B(3.1)=O.ODO 

B(3,2)=O.ODO 

B (3 , 3) = (BB1*K1*TL*TL*TL/ (6 . DO*XD) ) 
B(3,4)=O.ODO 

B (3 , 5) = (BB1*K1*TL*TL*TL/ (6 . DO*XD) ) 

B(4,1)=O.ODO 

B(4,2)=O.ODO 

B(4,3)=O.ODO 

B(4,4)=1.0D0 

B(4,5)=O.ODO 

B(5,1)=O.ODO 

B(5,2)=0.0D0 

B (5 , 3) =1 . ODO+ (BB2*K1*TL*TL*TL/ (6 . DO*XD) ) 
B(5,4)=O.ODO 

B (5 , 5) = (BB2*K1*TL*TL*TL/ (6 . DO*XD) ) 

RETURN 

END 
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C  PROGRAM  THESIS4.F0R  (Time  Delay-lst  Order  Approx  in  delta  ,T1) 
C  BIFURCATION  ANALYSIS 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2,L.NR,NV,IZ,MASS,NDELTA,NPSI,NY, 

&  NRDOT.NVDOT 

DIMENSION  A(4,4) ,B(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) ,ALFR(4) , 

&  ALFI(4),BETA(4),WR(4),WI(4) 

C 

OPEN  (11,FILE=>BIF1.RES’) 

OPEN  (12,FILE=’BIF2.RES’) 

OPEN  (13,FILE=’BIF3.RES’) 

C 

C  Vehicle  Parameters : 


IZ 

=0.0 

L 

=528 

RHO 

=1.94 

G 

=32.2 

XG 

=0.0 

MASS 

=0 . 0088 

U 

=3.0 

YRDOT  =  0.00000 
YVDOT  =-0.00912 
YR  =+0 . 00456 
YV  =-0 . 01434 

YPSI  =  0.01400 
YY  =  0.02000 

YDELTA=  0.00278 
NRDOT  =-0.00115 
NVDOT  =  0.00000 
NR  =-0 . 00296 
NV  =-0.00460 

NPSI  =  0.01000 

NY  =  -0.00250 

NDELTA=-0. 00139 

WRITE  (*.1001) 

READ  (*,♦)  WNMIN.WNMAX.IWN 

WRITE  (*.1002) 

READ  (*.*)  XDMIN.XDMAX.IXD 

WRITE  (*.1003) 
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ZETA 


READ  (*,*) 

WRITE(*,1100) 

READ  (*,*)  TL 
TL=TL*U/L 

DVR  =(IZ-NRDOT)*(MASS-YVDOT)- 
&  (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 

AAl 1= ( (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12= ( ( IZ-NRDDT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI ) /DVR 
AA22=( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) ♦YV) /DVR 
AA13=( (IZ-NRDOT) * (YR-MASS) + 

&  (YRDOT-MASS*XG)*(NR-MASS*XG))/DVR 

AA23= ( (NVDOT-MASS+XG) * (YR-MASS) + 

&  (MASS-YVDOT) *(NR-MASS*XG)) /DVR 

AA14=( (IZ-NRDOT) *YY+(YRDOT-MASS*XG) *NY)/DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BBl  =( (IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) ♦NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

EPS  =l.D-5 

ILMAX=1500 

DO  1  1=1, IWN 

WRITE  (*,2001)  I, IWN 

WN=WNMIN+(I-1)* (WNMAX-WNMIN)/ (IWN-1) 

K1=-WN**2.D0/BN0M 

K2=- (ANOM+2 . DO*ZETA*WN) /BNOM 

DO  2  J=1,IXD 

XD=XDMIN+( J-1) * (XDMAX-XDMIN) / (IXD-1) 

CALL  LINEAR (K1 , K2 , XD , AAll , AA12 , AA13 , AA14 , 

&  AA21,AA22,AA23,AA24,BB1,BB2,A,B,TL) 

CALL  RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , lERR) 

DO  11  IJE=1,4 

WR ( IJE) =ALFR (IJE) /BETA (I JE) 

WI (IJE) =ALFI (IJE) /BETA (IJE) 

CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

IF  (J.GT.l)  GO  TO  10 

DEOSOO=DEOS 

XDOO  =XD 


GO  TO  2 

DEOSNN=DEOS 

XDNN  =XD 

PR=DEOSNN*DEOSOO 

IF  (PR.GT.O.DO)  GO  TO  3 

LL=LL+1 

IF  (LL.GT.3)  STOP  1000 
IL=0 

XDO=XDOO 

XDN=XDNN 

DEOSO=DEOSOO 

DEOSN=DEOSNN 

XDL=XDO 

XDR=XDN 

DEOSL=DEOSO 

DEOSR=DEOSN 

XD=(XDL+XDR)/2.D0 

CALL  LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA2 1 , AA22 , 
AA23.AA24,BB1,BB2,A,B,TL) 

CALL  RGG(4,4,A,B,ALFR,ALFI,BETA,0,ZZZ,IERR) 

DO  12  IJE=1,4 

WR ( IJE) =ALFR ( I JE) /BETA ( I JE) 

WI ( I JE) =ALFI (IJE) /BETA ( I JE) 

CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS (XDL-XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 

GO  TO  4 

IF  (PRR.GT.O.DO)  STOP  3200 

XDO=XDM 

XDN=XDR 

DEOSO=DEOSM 

DEDSN=DEOSR 


IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 

4  LLL=10+LL 

WRITE  (LLL,*)  XD,WN 
3  XDOO=XDNN 

DEOSOO=DEOSNN 
2  CONTINUE 
1  CONTINUE 


C 


1001  FORMAT 

1002  FORMAT 

1003  FORMAT 
1100  FORMAT 
2001  FORMAT 

END 


(’  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN 
(»  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD 
(>  ENTER  DAMPING  RATIO') 

(’  ENTER  TIME  LAG  TL  (dimensional)’) 
(215) 


(dimensionless) ’) 
(dimensionless) ’) 


SUBROUTINE  DSTABL(DEOS,WR,WI, OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(4),WI(4) 

DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I) .LT.DEOS)  GO  TO  1 
DEOS=WR(I) 

IJ=I 

1  CONTINUE 
OMEGA=WI(IJ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 


SUBROUTINE  LINEAR(K1,K2,XD,AA11 ,AA12,AA13,AA14,AA21,AA22,AA23,AA24 
&  ,BB1,BB2,A,B,TL) 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DOUBLE  PRECISION  K1,K2 
DIMENSION  A(4,4),B(4,4) 

A(1,1)=O.ODO 

A(1,2)=O.ODO 

A(1,3)=1.0DO 

A(1,4)=O.ODO 

A(2,1)=AA11+BB1*K1-BB1*K1*TL/XD 
A (2 , 2) =AA12-BB1*K1*TL/XD 
A(2,3)=AA13+BB1*K2-BB1*K1*TL 
A(2 , 4) =AA14+BB1*K1/XD 
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A(3, 1)=AA21+BB2*K1-BB2*K1*TL/XD 
A (3 , 2) =AA22-BB2*K1 *TL/XD 
A (3 , 3) =AA23+BB2*K2-BB2*K1*TL 
A (3 , 4) =AA24+BB2*K1/XD 

A(4,1)=1.0D0 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

B(1,1)=1.0D0 

B(1,2)=0.0D0 

B(1,3)=0.0D0 

B(1,4)=0.0D0 

B(2,1)=0.0D0 

B(2,2)=1.0D0 

B(2,3)=BB1*K2*TL 

B(2,4)=0.0D0 

B(3,1)=0.0D0 
B(3,2)=0.0D0 
B (3 , 3) =1 . 0D0+ (BB2*K2*TL) 
B(3.4)=0.0D0 

B(4.1)=0.0D0 

B(4,2)=0.0D0 

B(4,3)=0.0D0 

B(4,4)=1.0D0 

RETURN 

END 
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C  PROGRAM  THESIS5.F0R  (Time  Delay-lst  Order  Approx  in  delta  &  y  ie.  in  T1  &  T2)) 
C  BIFURCATION  ANALYSIS 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2,L,NR,NV,IZ,MASS,NDELTA,NPSI,NY, 

&  NRDOT.NVDOT 

DIMENSION  A(4,4) ,B(4,4) ,FV1(4) .IV1(4) ,ZZZ(4,4) ,ALFR(4) , 

&  ALFI(4),BETA(4),WR(4) ,¥I(4) 

C 

OPEN  (11,FILE=’BIF1.RES') 

OPEN  (12,FILE=’BIF2.RESO 
OPEN  (13,FILE='BIF3.RES’) 

C 

C  Vehicle  Peirameters : 


IZ 

=0.0 

L 

=528 

RHO 

=1.94 

G 

=32,2 

XG 

=0.0 

MASS 

=0.0088 

U 

=3.0 

YRDOT 

=  0.00000 

YVDOT 

=-0,00912 

YR 

=+0.00456 

YV 

=-0.01434 

YPSI 

=  0.01400 

YY 

=  0.02000 

YDELTA=  0.00278 
NRDOT  =-0.00115 
NVDOT  =  0.00000 
NR  =-0.00296 
NV  =-0 . 00460 

NPSI  =  0.01000 

NY  =  -0 . 00250 
NDELTA=  -0.00139 

WRITE  (*,1001) 

READ  (*,*)  WNMIN,WNMAX,IWN 
WRITE  (*,1002) 

READ  ( * , *)  XDMIN , XDMAX , IXD 
WRITE  (*,1003) 

READ  (*,*)  ZETA 
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WRITE(*,1100) 

READ  (*,*)  TLl 

TL1=TL1*D/L 

WRITE(*,1101) 

READ  (*,*)  TL2 

TL2=TL2*U/L 

DVR  = ( IZ-NRDOT) * (MASS-YVDOT) - 
&  (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

AA11=( (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) +NPSI) /DVR 
AA12=( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13= ( (IZ-NRDOT) * (YR-MASS) + 

&  (YRDOT-MASS*XG)*(NR-MASS*XG))/DVR 

AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

&  (MASS-YVDOT) * (NR-MASS*XG) ) /DVR 

AA14=( (IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24= ( (NVDOT-MASS+XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BB 1  = ( ( IZ-NRDOT) *YDELTA- (MASS^XG-YRDOT) *NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS  =l.D-5 
ILMAX=1500 

DO  1  1=1, IWN 

WRITE  (*,2001)  I, IWN 

WN=WNMIN+ (I-l) * (WNMAX-WNMIN) / (IWN-1) 

Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- ( ANOM+2 . DO*ZETA*WN) /BNOM 
DO  2  J=1,IXD 

XD=XDMIN+ ( J- 1) * (XDMAX-XDMIN) / ( IXD- 1 ) 

CALL  LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , 

&  AA21,AA22,AA23,AA24,BB1,BB2,A,B,TL1,TL2) 

CALL  RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , lERR) 

DO  11  IJE=1,4 

WR(IJE)=ALFR(IJE) /BETA(IJE) 

WI (I JE) =ALFI (I JE) /BETA (IJE) 

CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 


IF  (J.GT.l)  GO  TO  10 

DEOSOO=DEOS 

XDOO  =XD 

LL=0 

GO  TO  2 

10  DEOSNN^DEOS 
XDNN  =XD 
PR=DEOSNN*DEOSOO 
IF  (PR.GT.O.DO)  GO  TO  3 
LL=LL+1 

IF  (LL.GT.S)  STOP  1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO^DEOSOO 
DEOSN=DEOSNN 
6  XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22, 
&  AA23 , AA24 , BBl , BB2 , A , B , TLl , TL2) 

CALL  RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , lERR) 

DO  12  IJE=1,4 

WR ( IJE) =ALFR ( I JE) /BETA ( IJE) 

WI (IJE) =ALFI (IJE) /BETA (IJE) 

12  CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDL''XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 

GO  TO  4 

5  IF  (PRR.GT.O.DO)  STOP  3200 
XDO=XDM 


67 


XDN=XDR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 

4  LLL=10+LL 

WRITE  (LLL,f)  XD,WN 
3  XDOO=XDNN 

DEOSOO=DEOSNN 
2  CONTINUE 
1  CONTINUE 


1001  FORMAT 

1002  FORMAT 

1003  FORMAT 

1100  FORMAT 

1101  FORMAT 
2001  FORMAT 

END 


(’  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN  (dimensionless)’) 
(’  ENTER  MIN,  MAX.  AND  INCREMENTS  OF  XD  (dimensionless)’) 
(’  ENTER  DAMPING  RATIO’) 

(’  ENTER  TIME  LAG  TLl  (dimensional)’) 

(’  ENTER  TIME  LAG  TL2  (dimensional)’) 

(215) 


C 


SUBROUTINE  DSTABL(DEOS,WR,WI, OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(4).WI(4) 

DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I) .LT.DEOS)  GO  TO  1 
DEOS=WR(I) 

IJ=I 

1  CONTINUE 
OMEGA=WI(IJ) 

OMEGA=DABS ( OMEGA) 

RETURN 

END 


SUBROUTINE  LINEAR(K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA2 1 , AA22 , AA23 , AA24 
&  ,BB1,BB2,A,B,TL1,TL2) 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DOUBLE  PRECISION  K1,K2 
DIMENSION  A(4,4) ,B(4,4) 

A(1,1)=O.ODO 

A(1,2)=0.0D0 

A(1,3)=1.0DO 

A(1,4)=0.0D0 
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A(2,1)=AA11+BB1*K1-BB1*K1*TL2/XD-BB1*K1*TL1/XD 
A(2 , 2) =AA12-BB1*K1*TL2/XD-BB1*K1*TL1/XD 
A(2,3)=AA13+BB1*K2-BB1*K1*TL1+BB1*K1*TL1*TL2/XD 
A(2,4)=AA14+BB1*K1/XD 

A(3 . 1)=AA21+BB2*K1-BB2*K1*TL2/XD-BB2*K1*TL1/XD 
A(3,2)=AA22-BB2*K1*TL2/XD-BB2*K1*TL1/XD 
A(3,3)=AA23+BB2*K2-BB2*K1*TL1+BB2*K1*TL1*TL2/XD 
A (3 , 4) =AA24+BB2*K1/XD 

A(4,1)=1.0D0 

A(4,2)=1.0D0 

A(4.3)=0.0D0 

A(4,4)=0.0D0 

B(1,1)=1.0D0 

B(1,2)=0.0D0 

B(1,3)=0.0D0 

B(1,4)=0.0D0 

B(2,1)=O.ODO 

B(2,2)=1.0D0-(BB1*K1*TL1*TL2/XD) 

B(2,3)=BB1*K2*TL1 

B(2,4)=0.0D0 

B(3,1)=0.0D0 

B(3,2)=-BB2*K1*TL1*TL2/XD 

B(3,3)=1.0D0+(BB2*K2*TL1) 

B(3,4)=0.0D0 

B(4,1)=0.0D0 

B(4,2)=O.ODO 

B(4,3)=O.ODO 

B(4,4)=1.0DO 

RETURN 

END 
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PROGRAM  HOPF.FOR 


Hopf  Bifurcation  Analysis 


Tliird  Order  Expansions:  First  Order  Approximation 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1.K2.K3,L,NR,NV,NDELTA,HPSI,NY,IZ,MASS, 

&  NRDOT , NVDOT , KIP , K2P , NVW , NVRR , NRVV , NYYY , NPPP 

DOUBLE  PRECISION  M11,M12,M13,M14,M21,M22,M23,M24, 

1  M31,M32,M33,M34,M41,M42.M43,M44, 

2  N11.N12,N13,N14,N21,N22,N23,N24, 

3  N3 1 , N32 , N33 , N34 , N41 , N42 , N43 , N44 , 

4  L21.L22,L23,L24,L31,L32,L33,L34, 

5  L41,L42,L43,L44 

DIMENSION  A (4, 4) ,T(4.4) ,TINV(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) 
DIMENSION  WR(4) ,W1(4) ,TSAVE(4,4) ,TLUD(4,4) ,IVLUD(4) ,SVLUD(4) 
DIMENSION  ASAVE(4,4) 


OPEN  (ll.FILE=’BIFl.RES’) 
OPEN  ( 15, FILE=’ HOPF. RES’) 


Vehicle  Pairameters : 


IZ 

=0.0 

L 

=528 

RHO 

=1.94 

G 

=32.2 

XG 

=0.0 

MASS 

=0.0088 

U 

=1.0 

YRDOT 

=  0.00000 

YVDOT 

=-0.00912 

YR 

=+0 . 00456 

YV 

=-0.01434 

YVW 

=-0.15391 

YVRR 

=-0.05476 

YRW 

=  0.04608 

YYYY 

=  0.46800 

YPPP 

=  0.00000 

YPSI 

=  0.01400 

YY 

=  0.02000 

YDELTA=+0 . 00278 


NRDOT  =-0.00115 
NVDDT  =  0.00000 
NR  =-0.00296 
NV  =-0.00460 
NVW  =-0.00336 
NVRR  =  0.00759 
NRW  =-0.05160 
NYYY  =  0.00000 
NPPP  =  0.00000 
NPSI  =  0.01000 
NY  =-0.00250 
NDELTA=-0.00139 


WRITE  (*,1003) 

READ  (*,*)  ZETA 

WRITE (*,1100) 

READ  (*,*)  TL 

TL=TL*U/L 
WRITE  (*,1006) 

READ  (*,*)  DO 


C 

DVR  =(IZ-NRDOT)*(MASS-YVDOT)- 
&  (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

AA11=( (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12= ( (IZ-NRDOT) *YV- (MASS*XG-yRDOT) *NV) /DVR 
AA21=( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13= ( ( IZ-NRDOT) * ( YR-MASS) + 

&  (YRDOT-MASS*XG)*(NR-MASS*XG))/DVR 

AA23=( (NVDOT-MASS*XG) * (YR-MASS) + 

&  (MASS-YVDOT) *(NR-MASS*XG)) /DVR 

AA14= ( ( IZ-NRDOT) *YY+ ( YRDOT-MASS*XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BBl  = ( ( IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 
C 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS  =l.D-5 
ILMAX=1500 
C 

IWN=1000 
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DO  1  11=1, IWN 

WRITE  (*,2001)  II 

READdl,*)  XD,WM 

Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- (ANOM+2 . DO*ZETA*WN) /BNOM 

K3=K2 


Start  Hopf  Bifurcation  Analysis 


Evaluate  Nonlinear  Rudder  Expansion  Coefficients 
K2=0.0D0 

K1P=K1-K1*TL*U/XD 

K2P=K2-K1*TL/XD 

DPPV=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1P*K2P 

+  0 . 5D0*K1*TL/XD  +  3.D0*K1*(TL*U)**3/(3.D0*XD**3) 
DPVV=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K2P*K2P 
+  3.D0*U*TL*TL*TL*K1/(3.D0*XD**3) 

DPPR=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1P*K3 
DPRR=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K3*K3 
DPPY=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1P*K1/XD 

-  3.D0*TL*TL*U*U*K1/(3.D0*XD**3) 

DPYY=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1*K1/ (XD**2) 

+  3.D0*TL*O*Kl/(3.DO*XD**3) 

DWR=-  ( 1 .  DO/  (3 .  DO*DO*  *2)  )  *3 .  D0*K2P*K2P*K3 
DVRR=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K2P*K3*K3 
DWY=-  ( 1 .  DO/  (3 . DO*DO**2)  )  *3 .  D0*K1*K2P*K2P/XD 

-  3.D0*K1*TL*TL/(3.D0*XD**3) 

DVYY=-  (1 . DO/ (3 . DO*DO**2) ) *3 . D0*K1*K1*K2P/ (XD**2) 

+  3.D0*TL*K1/(3.D0*XD**3) 

DRRY=-  (1 .D0/(3.D0*D0**2))*3.D0*K1*K3*K3/XD 
DRYY=-  (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1*K1*K3/ (XD**2) 

DPVR=-  ( 1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1P*K2P*K3 
DPVY=-  ( 1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1P*K1 *K2P /XD 

-  6.D0*TL*TL*U*K1/(3.D0*XD**3) 

DPRY=-  (1 . DO/ (3 . DO*DO**2) ) *6 . D0*K1P*K1*K3/XD 
DVRY=-  (1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1*K2P*K3 /XD 
DPPP=-  (1 . DO/ (3 . D0*D0**2) ) *1 . D0*K1P*K1P*K1P 

+  K1*TL*U/(6.D0*XD)  +  (K1*(TL*U)**3)/(3.D0*XD**3) 
DVW=-  (1 .  DO/  (3  .D0*D0**2)  )  *1 .  D0*K2P*K2P*K2P 
+  K1*(TL**3)/(3.D0*XD**3) 

DRRR=-  (1 . DO/ (3 . D0*D0**2) ) *1 . D0*K3*K3*K3 

DYYY=-  (1 . DO/ (3 .D0*D0**2) ) * (Kl/XD) **3-Kl/ (3 .D0*XD**3) 

-  K1/(3.D0*XD**3) 


Evaluate  Transformation  Matrix  of  Eigenvectors 


CALL  LINEAR(K1,K3.XD,AA11,AA12,AA13,AA14.AA21,AA22,AA23, 
&  AA24,BB1,BB2,A,TL) 

DO  11  1=1,4 
DO  12  J=l,4 

ASAVE(I,J)=A(I,J) 

12  CONTINUE 
11  CONTINUE 

CALL  RG(4.4,A,WR,WI,1,ZZZ,IV1,FV1,IERR) 

CALL  DSOMEG(IEV,WR,WI, OMEGA, CHECK) 

OMEGAO=OMEGA 
DO  50  1=1,4 

T(I,1)=  ZZZ(I,IEV) 

T(I,2)=-ZZZ(I,IEV+1) 

50  CONTINUE 

IF  (lEV.EQ.l)  GO  TO  13 

IF  (IEV.EQ.2)  GO  TO  14 

IF  (IEV.EQ.3)  GO  TO  15 

STOP  3004 

14  DO  60  1=1,4 

T(I,3)=ZZZ(I,1) 

T(I,4)=ZZZ(I,4) 

60  CONTINUE 
GO  TO  17 

15  DO  70  1=1,4 

T(I,3)=ZZZ(I,1) 

T(I,4)=ZZZ(I,2) 

70  CONTINUE 
GO  TO  17 

13  DO  16  1=1,4 

T(I,3)=ZZZ(I,3) 

T(I,4)=ZZZ(I,4) 

16  CONTINUE 

17  CONTINUE 
C 

C  Normalization  of  the  Critical  Eigenvector 

C 

CALL  NORMAL (T) 

C 

C  Definition  of  Mij 

C 

M11=T(1,1) 

M21=T(2,1) 

M31=T(3,1) 

M41=T(4,1) 

M12=T(1,2) 

M22=T(2,2) 
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M32=T(3,2) 

M42=T(4,2) 

M13=T(1,3) 

M23=T(2,3) 

M33=T(3,3) 

M43=T(4,3) 

M14=T(1.4) 

M24=T(2,4) 

M34=T(3,4) 

M44=T(4,4) 

C 

C  Definition  of  Lij 

C 

L21=  DPPV*M11*M11*M21  +  DPVV*M11*M21*M21  +  DPPR*M11*M11*M31 
&  +  DPRR*M11*M31*M31  +  DPPY*M11*M11*M41  +  DPYY*M11*M41*M41 

&  +  DVVR*M21*M21*M31  +  DVRR*M21*M31*M31  +  DVVY*M21*M21*M41 

k  +  DVYY*M21*M41*M41  +  DRRY*M31*M31*M41  +  DRYY*M31*M41*M41 

&  +  DPVR*M11*M21*M31  +  DPVY*M11*M21*M41  +  DPRY*M11*M41*M31 

k  +  DVRY*M21*M41*M31  +  DPPP*M11*M11*M11  +  DVVV*M21*M21*M21 

k  +  DRRR*M31*M31*M31  +  DYYY*M41*M41*M41 

C 

PPV  =  M11*M11*M22  +  2.0*M11*M12*M21 
PW  =  M12*M21*M21  +  2.0*M11*M21*M22 
PPR  =  M11*M11*M32  +  2.0*M11*M12*M31 
PRR  =  M12*M31*M31  +  2.0*M11*M31*M32 
PPY  =  M11*M11*M42  +  2.0*M11*M12*M41 
PYY  =  M41*M41*M12  +  2. 0*M11*M41*M42 
WR  =  M21*M21*M32  +  2.0*M31*M21*M22 
VRR  =  M22*M31*M31  +  2.0*M31*M32*M21 
VVY  =  M21*M21*M42  +  2.0*M41*M21*M22 
VYY  =  M22*M41*M41  +  2.0*M41*M42*M21 
RRY  =  M31*M31*M42  +  2.0*M41*M31*M32 
RYY  =  M32*M41*M4i  +  2.0*M41*M42*M31 
PVR  =  M11*M21*M32+M31*(M11*M22+M12*M21) 

PVY  =  M11*M21*M42+M41*(M11*M22+M12*M21) 

PRY  =  M11*M41*M32+M31*(M11*M42+M12*M41) 

VRY  =  M21*M41*M32+M31*(M21»M42+M22*M41) 

PPP  =  3.0*M11*M11*M12 
VW  =  3.0*M21*M21*M22 
RRR  =  3.0*M31*M31*M32 
YYY  =  3.0*M41*M41*M42 
C 

L22=DPPV*PPV+DPW=t=PW+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
k  +DWR*VVR+DVRR*VRR+DVVY*VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 

&  +DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVW*VVV 

&  +DRRR*RRR+DYYY*YYY 

C 

PPV  =  M12*M12*M21  +  2.0*M11*M12*M22 
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PVV  =  M11*M22*M22  +  2. 0*M12*M21*M22 
PPR  =  M12*M12*M31  +  2. 0*M11*M12*M32 
PER  =  M11*M32*M32  +  2 . 0*M12*M31*M32 
PPY  =  M12*M12*M41  +  2 . 0*M11*M12*M42 
PYY  =  M11*M42*M42  +  2. 0*M12*M41*M42 
WR  =  M22*M22*M31  +  2 . 0*M21*M22*M32 
VRR  =  M21*M32*M32  +  2 . 0*M22*M31*M32 
WY  =  M22*M22*M41  +  2 . 0*M21*M22*M42 
VYY  =  M21*M42*M42  +  2. 0*M22*M41*M42 
RRY  =  M32*M32*M41  +  2. 0*M31*M32*M42 
RYY  =  M31*M42*M42  +  2 . 0*M32*M41*M42 
PVR  =  M12*M22*M31  +  M32* (M11*M22+M12*M21) 

PVY  =  M12*M22*M41  +  M42* (M11*M22+M12*M21) 

PRY  =  M12*M42*M31  +  M32*(M11*M42+M12*M41) 

VRY  =  M22*M42*M31  +  M32*(M21*M42+M22*M41) 

PPP  =  3.0*M11*M12*M12 
WV  =  3.0*M21*M22*M22 
RRR  =  3.0*M31*M32*M32 
YYY  =  3.0*M41*M42*M42 

L23=DPPV*PPV+DPVV*PW+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
&  +DWR*WR+DVRR*VRR+DVVY*WY+DVYY*VYY+DRRY*RRY+DRYY*RYY 

&  +DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVW*WV 

&  +DRRR*RRR+DYYY*YYY 

L24=  DPPV*M12*M12*M22  +  DPVV*M12*M22*M22  +  DPPR*M12*M12*M32 
&  +  DPRR*M12*M32*M32  +  DPPY*M12*M12*M42  +  DPYY*M12*M42*M42 

&  +  DVVR*M22*M22*M32  +  DVRR*M22*M32*M32  +  DVVY*M22*M22*M42 

&  +  DVYY*M22*M42*M42  +  DRRY*M32*M32*M42  +  DRYY*M32*M42*M42 

&  +  DPVR*M12*M22*M32  +  DPVY*M12*M22*M42  +  DPRY*M12*M42*M32 

&  +  DVRY*M22*M42*M32  +  DPPP*M12*M12*M12  +  DVVV*M22*M22*M22 

&  +  DRRR*M32*M32*M32  +  DYYY*M42*M42*M42 

L31=L21 

L32=L22 

L33=L23 

L34=L24 

L21-L21*BB1*U*U 

L22=L22*BB1*U*U 

L23=L23*BB1*U*U 

L24=L24*BB1*U*U 

L31=L31*BB2*U*U 

L32=L32*BB2*U*U 

L33=L33*BB2*U*U 

L34=L34*BB2*U*U 

L41=-0 . (M21+U+M11/3 . 0) 

L42=-M11* (M12*M21+0 . 5*Mll*M22+0 . 5*U*M12*M11) 
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L43=-M12*(M11*M22+0.5*M12*M21+0.5*U*M11*M12) 

L44=-0 . 5*M12*M12* (M22+U*M12/3 . 0) 

C 

Cll=  (1/6 . 0)  *YVW*  (M21*M21*M21)+0 . 5*YVRR*  (M31*M31*M21)  + 

&  0.5*YBW*M21*M21*M31+(1/6.0)*YPPP*M11*M11*M11+(1/6.0)*YYYY*M41*M41*M41 

C12=(1/6.0)*YVW*(3*M21*M21*M22)+0.5*YVRR*(M31*M31* 

&  M22+2*M31*M32*M21)+0.5*YRVV*(M21*M21*M32+2*M31*M21*M22)+(1/6.0)*YPPP*3*M11*M11*M12 

&  +(1/6.0)*YYYY*3*M41*M41*M42 


C13=  ( 1/6 . 0)  *YVW*  (3*M21*M22*M22)  +0 . 5*YVRR*  (M21*M32* 

&  M32+2*M31*M32*M22)+0.5*YRVV*(M22*M22*M31+2*M32*M21*M22)+(1/6.0)*YPPP*3*M11*M12*M12 

&  +(1/6.0)*YYYY*3*M41*M42*M42 

C14=  ( 1/6 . 0)  *YVW*  (M22*M22*M22)  +0 . 5*YVRR*  (M32*M32*M22)  + 

&  0 . 5*YRW*M22*M22*M32+ (1/6.0) *YPPP*M12*M12*M12+ ( 1/6 . 0) *YYYY*M42*M42*M42 

C21=  ( 1/6 . 0)  *NVW*  (M21*M21*M21)+0 . 5*  (NVRR*M31*M31*M21)  + 

&  0.5*NRVV*M21*M21*M31+(l/6.0)*NPPP*Mll*Mll*Mll+(l/6.0)=t'NYYY*M41*M41*M41 

C22=(1/6.0)*NVW*(3*M21*M21*M22)+0.5*NVRR*(M31*M31* 

&  M22+2*M31*M32*M21)+0.5*NRW*(M21*M21*M32+2*M31*M21*M22)+(1/6.0)*NPPP*3*M11*M11*M12 

&  +(1/6.0)*NYYY*3*M41*M41*M42 

C23=  ( 1/6 . 0)  *NVW*  (3*M21*M22*M22)  +0 . 5*NVRR*  (M21+M32* 

&  M32+2*M31*M32*M22)+0 . 5*NRW* (M22*M22*M31+2*M32*M21*M22) + (1/6 . 0) ♦NPPP*3*M11*M12*M12 

&  +(1/6.0)*NYYY*3*M41*M42*M42 

C24=  (1/6.0)  *NVW*  (M22*M22*M22)  +0 . 5=t‘NVRR*M32*M32*M22+ 

&  0.5*NRW*M22*M22*M32+(1/6.0)*NPPP*M12*M12*M12+(1/6.0)*NYYY*M42*M42*M42 


D1 1= (Cll* (IZ-NRDOT) +C21* (YRDOT-M*XG) ) /DVR 
D12= (C12* (IZ-NRDOT) +C22* (YRDOT-M+XG) ) /DVR 
D13= (C13* (IZ-NRDOT) +C23* (YRDOT-M*XG) ) /DVR 
D14= (C14* (IZ-NRDOT) +C24* (YRDOT-M*XG) ) /DVR 

D21=(C11* (NVD0T-M*XG)+C21* (M-YVDOT) ) /DVR 
D22= (C12* (NVDOT-M*XG) +C22* (M-YVDOT) ) /DVR 
D23= (C13* (NVDOT-M*XG) +C23* (M-YVDOT) ) /DVR 
D24= (C14* (NVDOT-M*XG) +C24* (M-YVDOT) ) /DVR 
C 

C  Invert  Transformation  Matrix 

C 

DO  20  1=1,4 
DO  30  J=l,4 
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TINV(I,J)=0.0 

TSAVE(I,J)=T(I,J) 

CONTINUE 

CONTINUE 

CALL  DLUD(4,4,TSAVE,4,TLUD,IVLUD) 

DO  40  1=1,4 

IF  (IVLUD(I) .EQ.O)  STOP  3003 
CONTINUE 

CALL  DILU(4,4,TLUD,IVLUD,SVLUD) 

DO  80  1=1,4 
DO  90  J=l,4 

TINV(I,J)=TLUD(I,J) 

CONTINUE 

CONTINUE 

Check  Iiiv(T)*A*T 

IMULT=2 

IF  (IMULT.EQ.l)  CALL  MULT(TINV,ASAVE,T) 

IF  (IMULT.EQ.O)  STOP 

Definition  of  Nij 

Nll=TINV(i,l) 

N21=TINV(2,1) 

N31=TINV(3,1) 

N41=TINV(4,1) 

N12=TINV(1,2) 

N22=TINV(2,2) 

N32=TINV(3,2) 

N42=TINV(4,2) 

N13=TINV(1,3) 

N23=TINV(2,3) 

N33=TINV(3,3) 

N43=TINV(4,3) 

N14=TINV(1,4) 

N24=TINV(2,4) 

N34=TINV(3,4) 

N44=TINV(4,4) 

R11=N12*L21+N13*L31+N14*L41+N12*D11+N13*D21 

R12=N12*L22+N13*L32+N14*L42+N12*D12+N13*D22 

R13=N12*L23+N13=t‘L33+N14*L43+N12*D13+N13*D23 

R14=N12*L24+N13*L34+N14*L44+N12*D14+N13*D24 

R21=N22*L21+N23*L31+N24*L41+N22*D11+N23*D21 

R22=N22*L22+N23*L32+N24*L42+N22*D12+N23*D22 

R23=N22*L23+N23*L33+N24*L43+N22*D13+N23*D23 

R24=N22*L24+N23*L34+N24*L44+N22*D14+N23*D24 


c 

C  Evaluate  Alpha’  and  Omega’ 

C 

DINC=0.001 
XDR  =XD+DINC 
XDL  =XD-DINC 
XD  =XDR 

CALL  LINEAR (K1 , K3 , XD , AAl 1 . AA12 , AA13 , AA14 , AA21 , AA22 , AA23 , 
&  AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEOS 

OMEGR=FREQ 

XD=XDL 

C 

CALL  LINEAR (K1 , K3 , XD , AAl 1 , AA12 , AA13 , AA14 , AA21 , AA22 , AA23 , 
&  AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL (DECS , WR , WI , FREQ ) 

ALPHL=DEOS 

OMEGL=FREq 

C 

DALPHA= (ALPHR-ALPHL) / (XDR-XDL) 

DOMEGA= (OMEGR-OMEGL) / (XDR-XDL) 

C 

C  Evaluation  of  Hopf  Bifurcation  Coefficients 

C 

C0EF1=3 . 0*Rll+R13+R22+3 . 0*R24 
C0EF2=3 . 0*R21+R23-R12-3 . 0*R14 
AMU2  =-C0EFl/(8.0*DALPHA) 

BETA2=0.25*C0EF1 

TAU2  =- (C0EF2-D0MEGA*C0EF1/DALPHA) / (8 . O*0MEGA0) 

PER  =2.0*3. 1415927/OMEGAO 
PER  =PER*U/L 
C 

WRITE  (15,2002)  XD,WN,C0EF1 .DALPHA.PER 


1  CONTINUE 
C 

1001  FORMAT  (’  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN  (dimensionless)’) 

1002  FORMAT  (’  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD  (dimensionless)’) 

1003  FORMAT  (’  ENTER  DAMPING  RATIO’) 
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1006 

FORMAT 

(’  ENTER 

DSAT 

0 

1100 

FORMAT 

(>  ENTER 

TIME 

LAG 

2001 

FORMAT 

(215) 

2002 

FORMAT 

(5E15.5) 

END 

SUBROUTINE  DSTABL(DEOS,WR,WI, OMEGA) 

C 

C  EVALUATES  THE  EIGENVALUE  WITH  THE  MAXIMUM  REAL  PART 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  WR(4),WI(4) 

DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I) .LT.DEOS)  GO  TO  1 
DEOS=WR(I) 

IJ=I 

1  CONTINUE 
OMEGA=WI(IJ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 

C 

SUBROUTINE  LINEAR (K1,K2,XD,AA11,AA12,AA13,AA14, AA2 1 , 

&  AA22,AA23.AA24,BB1.BB2.A.TL) 

C 

C  FORMS  THE  LINEARIZED  MATRIX  A  (time  delay  1st  order  approximation 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DOUBLE  PRECISION  K1,K2 
DIMENSION  A (4, 4) 

A(1,1)=O.ODO 

A(1,2)=O.ODO 

A(1,3)=1.0DO 

A(i,4)=O.ODO 

A(2,1)=AA11+BB1*K1-BB1*K1*TL/XD 

A(2 , 2)  =AA12-BB1=»K1*TL/XD 

A(2,3)=AA13+BB1*K2 

A (2 , 4) =AA14+BB1*K1/XD 

A (3 , 1) =AA21+BB2*K1-BB2*K1*TL/XD 

A (3 , 2) =AA22-BB2*K1*TL/XD 

A(3,3)=AA23+BB2*K2 

A (3 . 4) =AA24+BB2*K1/XD 

A(4,1)=1.0D0 

A(4,2)=1.0DO 
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A(4,3)=0.0D0 

A(4.4)=0.0D0 

RETURN 

END 


C. 


SUBROUTINE  DSOMEG(IJK.WR,WI, OMEGA , CHECK ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(4),WI(4) 

CHECK=-1.0D+25 
DO  1  1=1,4 

IF  (WR(I) .LT. CHECK)  GO  TO  1 
CHECK=WR(I) 

IJ=I 

1  CONTINUE 

OMEGA=DABS(WI(IJ)) 

IF  (WI(IJ).GT.O.DO)  IJK=IJ 
IF  (WI(IJ).LT.O.DO)  IJK=IJ-1 
RETURN 
END 


SUBROUTINE  NORMAL (T) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  T(4,4) ,TNOR(4,4) 

CFAC=T (1,1) **2+T (1 , 2) **2 

IF  (DABS(CFAC) .LE. (l.D-10))  STOP  4001 

TN0R(1,1)=1.D0 

TN0R(2,1)=(T(1,1)*T(2,1)+T(2,2)*T(1,2))/CFAC 

TN0R(3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2))/CFAC 

TN0R(4.1)=(T(1,1)*T(4,1)+T(4,2)*T(1,2))/CFAC 

TN0R(1,2)=0.D0 

TN0R(2,2)=(T(2,2)*T(1.1)-T(2,1)*T(1,2))/CFAC 
TN0R(3,2)=(T(3,2)*T(1,1)-T(3,1)*T(1,2))/CFAC 
TN0R(4,2)=(T(4,2)*T(1,1)-T(4,1)*T(1,2))/CFAC 
DO  1  1=1,4 
DO  2  J=l,2 

T(I,J)=TNOR(I,J) 

2  CONTINUE 
1  CONTINUE 
RETURN 
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END 


SUBROUTINE  MULT(TINV,A,T) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  TINV(4,4) ,A(4,4) ,T(4,4) ,A1(4,4) ,A2(4,4) 
DO  1  1=1,4 
DO  2  J=l,4 
Aid,  J)=O.DO 
A2(I,J)=0.D0 

2  CONTINUE 
1  CONTINUE 

DO  3  1=1,4 
DO  4  J=l,4 
DO  5  K=l,4 

A1(I,J)=A(I,K)*T(K,J)+A1(I,J) 

5  CONTINUE 
4  CONTINUE 

3  CONTINUE 
DO  6  1=1,4 

DO  7  J=l,4 
DO  8  K=l,4 

A2(I,J)=TINV(I,K)*A1(K,J)+A2(I,J) 

8  CONTINUE 

7  CONTINUE 

6  CONTINUE 

DO  11  1=1,4 

C  WRITE  (*,101)  (A(I,J) ,J=1,4) 

11  CONTINUE 
DO  12  1=1,4 

C  WRITE  (*,101)  (T(I, J) ,J=1,4) 

12  CONTINUE 
DO  10  1=1,4 

WRITE  (*,101)  (A2(I,J) ,J=1,4) 

10  CONTINUE 

C  WRITE  (*,101)  A2(l,l) 

RETURN 

101  FORMAT  (4E15.5) 

END 


81 


LIST  OF  REFERENCES 


1.  Beck,  R.  F.  [1976]  “Forces  and  moments  on  a  ship  moving  in  a  canal”, 
Report  No.  179,  Dept,  of  Naval  Architecture  and  Marine  Engineering, 
The  University  of  Michigan,  Ann  Arbor. 

2.  Bernitsas,  M.  M.  and  Kekridis,  N.  S.  [1984]  “Nonlinear  simulation  of  time 
dependent  towing  of  ocean  vehicles” ,  Report  No.  283,  Dept,  of  Naval 
Architecture  and  Marine  Engineering,  The  University  of  Michigan,  Ann 
Arbor. 

3.  Chow,  S.-N.  and  Mallet-Paret,  J.  [1977]  “Integral  averaging  and  bifur¬ 
cation”,  Journal  of  Differential  Equations,  26,  pp.  112-159. 

4.  Clayton,  B.  R.  and  Bishop,  R.  E.  D.  [1982]  Mechanics  of  Marine  Vehicles 
(Gulf  Publishing  Company,  Houston). 

5.  Comstock,  J.  P.  [1977]  Principles  of  Naval  Architecture  (The  Society  of 
Naval  Architects  and  Marine  Engineers,  New  York). 

6.  Guckenheimer,  J.  and  Holmes,  P.  [1983]  Nonlinear  Oscillations,  Dynam¬ 
ical  Systems,  and  Bifurcations  of  Vector  Fields  (Springer- Verlag,  New 
York). 

7.  Hassard,  B.  and  Wan,  Y.H.  [1978]  “Bifurcation  formulae  derived  from 
center  manifold  theory” ,  Journal  of  Mathematical  Analysis  and  Applica¬ 
tions,  63,  pp.  297-312. 

8.  Venne,  D.  V.  [1992]  “Effects  of  positional  information  time  lags  on  motion 
stability  of  autonomous  vehicles”.  Master’s  Thesis,  Naval  Postgraduate 
School,  Monterey. 


83 


84 


INITIAL  DISTRIBUTION  LIST 


No.  Copies 

1.  Defense  Technical  Information  Center  . 2 

8725  John  J.  Kingman  Rd.  STE  0944 

Ft.  Belvoir,  VA  22060-6218 

2.  Dudley  Knox  Library  . 2 

Naval  Postgraduate  School 

411  Dyer  Rd. 

Monterey,  California  93943-5101 

3.  Chairman,  Code  ME . 1 

Department  of  Mechanical  Engineering 

Naval  Postgraduate  School 
Monterey,  California  93943 

4.  Professor  Fotis  A.  Papoulias,  Code  ME/Pa . 3 

Department  of  Mechanical  Engineering 

Naval  Postgraduate  School 
Monterey,  California  93943 

5.  Panos  E.  Kapasakis . 3 

Argyriou  Mathesi  4 

Salamina  18900 
Greece 

6.  Embassy  of  Greece . 1 

Naval  Attache 

2228  Massachusetts  Avenue,  N.W. 

Washington,  D.C.  20008 

7.  Naval  Engineering  Curricular  Office,  Code  34 . 1 

Naval  Postgraduate  School 

Monterey,  California  93943 


85 


